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Abstract.  This  thesis  concerns  the  use  of  time-split  methods  for  the  numerical  solu¬ 
tion  of  time-dependent  partial  differential  equations.  Frequently  the  differential  operator 
splits  additivcly  into  two  or  more  pieces  such  that  the  corresponding  subproblctns  are 
each  easier  to  solve  than  the  original  equation,  or  arc  best  handled  by  different  tech¬ 
niques.  In  the  time-split  method  the  solution  to  the  original  equation  is  advanced  by  al¬ 
ternately  solving  the  subproblems.  In  this  thesis  a  unified  approach  to  splitting  methods 
is  developed  which  simplifies  their  analysis.  Particular  emphasis  is  given  to  splittings  of 
hyperbolic  problems  into  subproblcms  with  disparate  wave  speeds. 

Three  main  aspects  of  the  method  are  considered.  The  first  is  the  accuracy  and 
efficiency  of  the  time-split  method  relative  to  unsplit  methods.  We  derive  a  general 
expression  for  the  splitting  error  and  use  it  to  compute  the  overall  truncation  error  for 
the  time-split  method.  This  is  then  used  to  analyze  its  efficiency,  measured  by  the  amount 
of  work  required  to  obtain  a  given  accuracy. 

^ The  second  topic  is  stability  for  split  methods.  After  a  demonstration  that  in 
general  the  product  of  two  stable  operators  need  not  be  stable,  some  important  classes  of 
hyperbolic  splittings  are  identified  for  which  the  product  of  stable  approximate  solution 
operators  is  in  fact  stable. 

<~'The  final  topic  is  the  proper  specification  of  boundary  data  for  the  intermediate  solu¬ 
tions,  c.g.,  the  solution  obtained  after  solving  only  one  of  the  subproblems.  A  procedure 
is  described  which,  for  many  problems,  can  be  used  to  transform  the  given  boundary 
conditions  for  the  original  equation  into  arbitrarily  accurate  boundary  coalitions  for  the 
intermediate  solutions.  Stability  of  the  initial-boundary  value  problem  Is  also  discussed. 
^  ^The  main  emphasis  is  on  hyperbolic  problems,  and  the  one-dimensional  shallow 
water  equations  are  used  as  a  specific  example  throughout.  The  final  chapter  is  devoted  to 
some  other  applications  of  the  theory.  Two-dimensionaf  hyperbolic  problems,  convection- 
diffusion  equations,  and  the  Peaceman-Rachford  ADI  method  for  the  heat  equation  are 
considered. 
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1.  Introduction 


1.1.  Ownrfaw. 

Splitting  methods  of  one  form  or  another  are  frequently  used  in  computing  numerical 
solutions  to  partial  differential  equations.  This  thesis  concerns  one  wide  class  of  splitting 
methods  which  will  be  referred  to  as  time-aplit  method a.  Such  methods  are  also  known  as 
fractional  ate p  methoda.  These  methods  apply  to  time-dependent  equations  of  the  form 
ut  =  X(u)  for  which  the  differential  operator  A  splits  additiveiy  into  two  or  more  pieces, 
say  A(u)  —  Ai{u)  +  y?a(u),  such  that  the  subproblems 

Ut  =  *i(u) 

and 

«t  =  A  a(t») 

are  each  easier  to  solve  than  the  original  problem,  or  are  best  handled  by  different  tech¬ 
niques.  In  the  time-split  method,  the  solution  to  the  original  problem  is  advanced  by 
alternating  between  (approximately)  solving  each  of  the  two  subproblcms.  For  example, 
a  multi-dimensional  problem  may  be  split  into  one-dimensional  subproblems,  convection- 
diffu»'»n  or  the  Navier-Stokes  equations  may  be  split  into  hyperbolic  and  parabolic  sub- 
prablcnib,  or  a  purely  hyperbolic  problem  may  be  split  into  subproblcms  with  disparate 
wave  speeds. 

The  aims  of  this  thesis  are  twofold.  The  first  is  to  present  a  unified  framework 
for  studying  various  aspects  of  time-split  methods.  The  main  idea  is  to  decompose  the 
derivation  of  a  time-split  method  into  two  steps.  First  the  exact  solution  operator  for 
the  original  problem  is  approximated  to  second  order  accuracy  by  a  product  of  exact 
solution  operators  for  the  subproblcms.  Then  these  exact  solution  operators  arc  replaced 
by  second  order  accurate  numerical  approximations.  Many  commoniy-uscd  splitting 
methods  can  be  viewed  in  this  manner  (see  Section  1.6).  This  viewpoint  is  not  new,,  but 
some  of  Its  consequences  have  not  been  fully  exploited. 

One  advantage  of  this  approach  is  that  the  errors  in  the  resulting  numerical  ap¬ 
proximation  can  be  decomposed  into  errors  due  to  splitting  the  exact  solution  operator 
and  errors  due  to  numerically  solving  the  subproblcms.  Ttic  latter  errors  arc  well  under¬ 
stood  when  standard  numerical  methods  are  applied.  Section  2.3  contains  some  general 
expressions  for  the  splitting  error.  This  decomposition  of  errors  aids  in  analysing  the 
efficiency  of  the  time-split  method,  defined  as  the  amount  of  work  required  to  obtain 
a  given  accuracy.  The  sise  of  the  splitting  error  relative  to  the  truncation  errors  of 
the  numerical  methods  employed  plays  a  critical  role  in  determining  the  optimal  choice 
of  mesh  ratio  for  the  time-split  method  and  in  determining  how  efficient  the  resulting 
method  will  be  relative  to  uusplit  methods.  This  is  investigated  in  detail  in  Chapter  2. 
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Another  advantage  of  this  viewpoint  is  that  the  intermediate  solutions  (e.g.,  the  solu¬ 
tion  obtained  after  solving  only  one  of  the  subproblcms)  take  on  physical  meaning.  They 
are  second  order  accurate  approximate  solutions  to  some  differential  equation  (though 
not  to  the  original  equation).  This  is  an  important  realisation,  particulary  when  we 
attempt  to  specify  boundary  conditions  for  the  intermediate  solutions.  Such  boundary 
conditions  are  often  required  to  implement  the  time- split  method  and  have  frequently 
been  specified  in  an  ad  hoc  manner,  e.g.,  the  boundary  conditions  from  the  original  equa¬ 
tion  are  imposed  on  the  intermediate  solutions  as  well.  More  sophisticated  approaches, 
such  as  the  method  of  undetermined  functions[56],  derive  correct  boundary  conditions 
based  on  the  finite  difference  equations.  However,  by  viewing  the  intermediate  solution 
as  an  approximation  to  a  differential  equation,  it  is  often  possible  to  derive  appropriate 
boundary  conditions  without  regard  to  the  finite  difference  methods  employed.  The  given 
boundary  conditions  for  the  original  equation  are  transformed  into  boundary  conditions 
appropriate  for  the  oubproblems.  This  is  the  subject  of  Chapter  4. 

The  second  aim  of  this  thesis  is  to  investigate  the  applicability  of  the  time-split 
method  to  one  particular  class  of  problems,  namely  to  hyperbolic  systems  of  equations 
which  are  split  into  subproblcms  with  disparate  wave  speeds.  The  original  problem 
either  has  all  fast  waves  or  some  fast  waves  and  some  slow  waves.  This  splitting  may  be 
advantageous  if  the  “fast”  subproblem  can  be  solved  more  efficiently  than  the  full  system. 
The  remaining  subproblem  can  also  be  solved  more  efficiently  than  the  full  system  since 
only  slow  waves  are  present.  Such  problems  are  described  in  detail  in  Section  1.4. 

Time-split  methods  for  hyperbolic  problems  have  not  been  studied  extensively  in 
the  past,  but  the  results  presented  here  indicate  tkat  in  many  situations  they  are  quite 
valuable. 

Hyperbolic  problems  also  provide  specific  examples  for  the  general  theory  being 
developed.  For  example,  both  the  efficiency  analysis  of  the  time-split  method  and  the 
procedure  for  specifying  intermediate  boundary  conditions  arc  introduced  by  considering 
hyperbolic  examples.  A  few  other  applications  are  treated  in  Chapter  5. 

1.2.  Some  partial  differential  equations  and  finite  difference  methods. 

Time-dependent  partial  differential  equations  arise  in  modeling  a  wide  variety  of 
physical  phenomena.  Simple  examples  in  two  space  dimensions  include  the  parabolic 
heat  equation 

=  u„  +  Uyy  (1.1) 

and  the  hyperbolic  wave  equation 

utt  =  u„  +  tivw.  (1.2) 

The  latter  equation  can  be  rewritten  as  a  first  order  hyperbolic  system  of  equations  in 
the  variables  ut,  ux,  and  uy.  A  general  first  order  hyperbolic  system  has  the  form 

t*i  —  Au,  +  Buy  (1.3) 

where  u  is  now  a  vector  and  A  and  B  are  diagonaliiablc  matrices  with  real  eigenvalues. 
For  the  wave  equation  (1.2)  A  and  B  arc  constant,  but  in  a  more  general  variable 

coefficient  problem,  A  and  B  could  depend  on  z,  y,  and  t.  If  A  and  B  also  depend 
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on  the  solution  «  then  the  problem  la  said  to  be  qua si/inear.  The  inviscid  Euler  equations 
of  gas  dynamics  are  of  this  form,  for  example. 

Practical  problems  often  include  both  first  and  second  order  spatial  derivatives. 
The  simplest  example  is  the  scalar  convection-  diffusion  equation  in  one  space  dimension, 
which  has  the  form 

Ut  me  eum  +  UUXX  (1.4) 

for  some  constants  e  and  v  >  0.  More  realistically,  the  compressible  Navier-Stokes 
equation t  for  viscous  flow  in  two  dimensions  constitute  a  quasilinear  system  of  the  form 

=  Aua  +  Buv  +  CuIX  +  Du,*  +  25u«,  (1.5) 

where  each  of  the  matrices  is  a  function  of  u. 

Inhomogeneous  terms  can  also  arise  in  practice.  For  example,  the  primitive  equations 
of  atmospheric  flow  (the  shallow  water  equations)  are  a  quasilinear  system  of  the  form 
(1.3)  with  an  undifferentiated  vector  F{u)  representing  Coriolis  forces  added  to  the  right 
hand  side. 

Lower-order  terms  also  occur  in  reaction- diffusion  equations  of  the  form 


ut  =  u(uxx  +  uvv)  +  F(u).  (1.6) 

Here  F  represents  the  chemical  kinetics  of  a  reacting  system  with  diffusivity  v  >  6. 

All  of  the  examples  given  above  are  of  the  general  form 

ut  =  A[u)  (1.7) 

where  A(u)  depends  on  u  and  its  spatial  derivatives.  It  may  also  depend  on  t  and  the 
spatial  variables  although  this  dependence  is  not  explicitly  shown. 

Initial  boundary  value  problems  will  be  discussed  in  detail  later  in  this  thesis,  but  for 
now  we  restrict  our  attention  to  the  Cauchy  problem,  which  consists  of  the  equation  (1.7) 
on  the  unbounded  spatial  domain  — oo  <  x  <  oo,  — oo  <  y  <  oo  (in  two  dimensions) 
together  with  initial  data  u(x,y,to)  =  f[x,y). 

If  the  problem  is  well-posed  (as  all  of  the  examples  above  are)  then  the  initial 
conditions  uniquely  determine  the  solution  at  any  later  time  t|.  We  write 


«(t|)  =  S{ti,to)u(to). 


(1.8) 


In  general  the  solution  operator  »9(ti,  to)  is  nonlinear,  but  satisfies  the  semigroup  property 

5(t8,t0)  =  5(<2,t,)5(t,fio)  (1.9) 

if  t0  <  <  t9. 

if  t  docs  not  appear  explicitly  in  the  coefficients  of  the  differential  equation,  then  the 
equation  is  said  to  be  autonomous  and  the  solution  operator  depends  only  on  the  time 
■  elapsed: 

S(t,,t0)  =  S(fi-to). 


For  notationat  convenience  we  will  assume  that  this  is  so  unless  otherwise  stated. 


Most  practical  problems  cannot  be  solved  exactly.  Instead  the  solution  must  be  ap¬ 
proximated  numerically.  We  will  be  concerned  only  with  finite  difference  approximations. 
For  such  methods  a  grid  is  laid  out  over  the  spatial  domain  and  an  approximate  solution 

at  all  gridpoints  is  obtained  at  each  of  a  sequence  of  times  to>  ti .  In  general  we 

assume  that  to  =  0  and  that  t„  =  nk  for  some  times tep  k.  For  convenience  we  assume 
that  the  grid  is  uniform  with  equal  mesh  spacing  h  in  all  spatial  coordinate  directions, 
although  this  is  not  necessary.  We  will  always  use  X  to  denote  the  mtth  ratio: 

X  =  k/h. 

Numerical  approximations  are  denoted  by  capital  letters.  In  one  space  dimension 
U *  is  the  approximation  to  u (xm,tn)  where  xm  =  mh.  In  higher  dimensions  more 
subscripts  are  added. 

We  will  restrict  our  attention  to  two-level  difference  schemes.  This  means  that  l/n+1 
is  determined  solely  by  Un  via  some  relation 

Un+i  =  Q(k)Un.  (1.10) 


This  is  the  difference  analogue  to 


«(t„+i)  =  S(k)u(tn) 

and  the  finite  difference  operator  Q(k)  is  an  approximation  to  the  solution  operator 
S(k).  The  method  is  said  to  be  accurate  of  order  p  if,  for  smooth  functions  it,  the  local 
truncation  error  ( Q(k )  —  S(k))u  is  0(kp+x)  as  *  0  with  some  fixed  relation  between  k 

and  h. 

As  an  example,  consider  a  one-dimensional  constant  coefficient  hyperbolic  equation 

Ut  =  Aux. 

Here  u  6  i»  a  vector  and  A  €  fftrXr  is  a  square  matrix.  By  Taylor  scries  expansions 
we  find  that  the  exact  solution  satisfies 

u(z,  t  +  k)  =  u(x,  t)  +  kut(x,  <)  +  j  k^uuix,  <)  +  ••• 

=  u(z,  t)  +  kAux(x,  t)  +  JfcaAsu«,(z,  <)  +  ••• 

=  (I  +  kAdx  +  lk2A*dl  +  "-)u{x,t) 

—  exp(kAdx)u(x,  <). 

We  thus  have  S(k)  =  r\p(kAOx),  as  defined  by  the  series  expansion  for  the  exponential. 
It  is  convenient  to  use  this  exponential  notation  for  the  solution  operators  of  constant 
coefficient  problems. 

The  Lax- Wen  dr  off  method.  If  the  expansion  for  exp  (kAdx)  is  truncated  after 
the  first  three  terms  and  the  differential  operators  dx  and  d\  are  replaced  by  appropriate 
finite  difference  operators,  we  obtain  the  familiar  Lax-Wendroff  method : 

fC+1  =(T  +  kAD0+  {k'A'D+DJph  (1.11) 
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where 


D°t/«  =  ~  tf—i), 

D+Um  =*  —  Um), 

D-Um=±(Um-Um- 1), 

D+D-Um  =  -  2l/m  +  t/m-l). 

The*  value  1/JJ+1  it  thus  determined  by  the  values  f/JJ,,  and  This  is 

conveniently  denoted  by  showing  the  stencil  of  the  the  method  as  in  Figure  1.1. 


FIG.  1.1.  Stencil  for  Lax-Wendroff. 


The  numerical  operator  Q(k)  is  defined  by  equation  (1.11).  This  Lax-Wendroff 
operator  appears  so  often  in  the  sequel  that  we  will  introduce  the  following  notation 
for  it,  which  shows  the  dependence  on  the  coefficient  matrix  A  explicitly: 

LW(A,  k)  =s  I  +  kAD0  +  \k*A%D+D_.  (1.12) 

Strictly  speaking,  this  operator  also  depends  on  h,  or,  equivalently,  on  the  mesh  ratio 
X  =  k/h,  but  X  will  be  assumed  to  be  fixed.  Analogous  methods  can  be  defined  for 
variable  coefficient  or  quasilincar  hyperbolic  systems.  The  same  generic  symbol  LW[A,  k) 
will  be  used  Tor  all  of  these  methods  although  in  general  they  will  be  more  complicated 
than  in  (1.12). 

The  Lax-Wendroff  method  is  second  order  accurate  since  the  local  truncation  error 
is  0(k3): 


\LW{A,  k)  -  exp(M0,)ju(x,  t)  =  -  J*S(AS  -  &A)«„S  +  0(k4).  (1.13) 


The  Crank- Nicolson  method.  As  another  example,  consider  the  one  dimensional 
heat  equation 

o,  =  (1.14) 

The  solution  operator  for  this  problem  is  S(k)  =  exp(/td®).  Explicit  methods  for 
parabolic  problems  are  generally  stable  only  if  the  Limcstep  k  is  very  small  relative  to 
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h.  For  this  reason  implicit  methods  are  often  used  instead.  One  popular  method  is  the 
second  order  accurate  Crank- Nicolton  method, 

(1  -  J*D+D_)l/"+1  =  (1  +  \kD+D-)U*mt  (1.15) 

for  which 

Q(k)  =  (1  -  bkD+D-)~l(  1  + 

This  corresponds  to  using  a  rational  approximation  to  the  exponential  solution  operator. 
To  implement  this  method  a  tridiagonal  system  of  linear  equations  must  be  solved  at  each 
iteration.  This  can  be  done  quite  efficiently.  Because  all  of  the  must  be  determined 
simultaneously,  the  method  is  said  to  be  implicit  Lax-Wendroff,  by  comparison,  is  an 
explicit  method.  The  stencil  for  Crank- Nicolson  is  shown  in  Figure  1.2. 


FIG.  1.2.  Stencil  For  Crank-Nicolson. 


In  two  space  dimensions  the  heat  equation  (1.1)  can  be  solved  by  a  similar  method: 

(1  -  **(D+,D_,  +  D+VD-,))UXJ  =  (1  +  ifc(D+x»-,  +  Z>+si>-„))tC,y,  (1.16) 

where,  for  example,  D+x  is  the  forward  difference  operator  in  the  z-direction.  Unfortun¬ 
ately,  this  no  longer  leads  to  a  tridiagonal  system  of  equations  but  rather  to  a  more  com¬ 
plicated  system  which  cannot  be  solved  nearly  as  efficiently.  It  was  this  problem  which 
led  to  the  introduction  of  some  of  the  first  splitting  methods.  One  such  method  is  the 
locally  one-dimensional  (LOD)  method  in  which  the  solution  to  (1.1)  is  advanced  by  first 
solving  u t  =  u„  approximately  using  (1.15)  and  then  solving  ut  =  uvy  approximately 
using  the  same  method  in  the  y-direction.  In  this  manner  only  one-dimensional  problems 
need  be  solved.  The  LOD  method  is  one  special  case  of  the  time-split  method,  which  will 
now  be  described  more  generally. 

1.3.  The  time-split  method. 

Consider  again  the  general  problem  (1.7)  and  suppose  that  the  function  >i(u)  splits 
additivcly  into  two  or  more  pieces  which  arc  most  naturally  handled  separately.  Restricting 
our  attention  to  two  pieces,  suppose  A  is  of  the  form 

A(u)  =  *i(u)  +  *a(u),  (1.17) 

where  each  or  the  subproblems 

ut  =  A  i(u)  (1.18a) 
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and 


(1.18b) 


is  easier  to  solve  than  the  full  problem  (1.7).  As  we  have  already  seen,  this  is  the  case  for 
the  heat  equation  (1.1)  when  ^i(u)  =  u„  and  ^2(ti)  =  uyy.  It  may  also  prove  useful  to 
handle  the  different  space  dimensions  in  the  hyperbolic  system  (1.3)  separately.  For  other 
equations  the  natural  splitting  is  between  terms  describing  different  physical  processes.  In 
the  convection-diffusion  equation  (1.4)  we  may  take  Ai(u)  =  cux  and  ^2(u)  =  t/u„,  thus 
splitting  the  mixed  problem  up  into  separate  hyperbolic  and  parabolic  equations.  The 
reaction-diffusion  equation  (1.6)  might  be  handled  similarly.  The  Navier-Stokes  equation 
(1.5)  could  well  be  split  into  five  separate  pieces. 

Splittings  have  long  been  used  for  all  of  these  problems  and  in  many  other  contexts 
as  well.  Some  history  and  references  are  given  in  Section  1.6. 

We  now  discuss  in  more  detail  the  implementation  of  the  time-split  method  once  a 
splitting  of  the  form  (1.17)  has  been  decided  upon.  The  subproblems  (1.18a)  and  (1.18b) 
have  corresponding  solution  operators  Si(k)  and  S2(fc).  The  basic  assumption  is  that 
these  operators  are  easier  to  approximate  than  S(k)  is.  The  time-split  method  is  based 
on  the  fact  that 

S{k)  «  S2(jfc)S,(*)  (1.19) 

when  k  is  small.  In  some  cases  this  splitting  is  in  fact  exact.  For  the  heat  equation  (1.1) 
with  the  LOD  splitting,  for  example,  we  have  S(k)  =  exp(ifc(3|  +  d*))  while  Si(ifc)  = 
exp(fci9®),  S2(k)  =  exp (kdy).  Since  the  differential  operators  d\  and  dy  commute,  we 
find  that  S(k)  =  S2(fc)Si(A:).  For  variable  coefficient  problems,  or  systems  or  equations, 
the  splitting  (1.19)  is  not  exact  in  general.  For  example,  the  same  LOD  splitting  on  the 
constant  coefficient  hyperbolic  system  (1.3)  has  an  error 


Si{k)S\{k)  —  S(k)  —  exp(fc#dy)cxp(fcvic)x)  —  exp^Aef*  +  Bdx)) 
=  ±k\BA  -  AB)dxdy  +  0{k3) 


(1.20) 


as  can  be  verified  by  expanding  the  exponentials.  In  this  case  the  splitting  is  exact  only  if 
the  matrices  A  and  B  commute.  Otherwise  the  local  error  on  smooth  solutions  is  0(k2) 
and  hence  the  splitting  is  only  first  order  accurate. 

A  simple  second  order  splitting  was  introduced  for  this  problem  by  Strang[49]  who 
noted  that 

exp(fc(/td*  +  Bdy))  =  exp(5<:A9I)exp(fcI7^y)exp(^fcAc)I)  +  0(k3). 

In  fact  the  same  type  of  splitting  is  second  order  accurate  (on  smooth  solutions)  for 
general  problems  of  the  form  (1.7).  The  general  Strang  splitting  is 

S(fc)«S,(ifc/2)S2()fc)Si(A/2).  (1.21) 

If  the  equation  depends  explicitly  on  t,  then  the  appropriate  form  of  the  splitting  is 

S(t  +  k,  t)  Si(t  +k,t+  kk)S-i(t  +  k,  t)St{t  +  Ajfc,  t). 

fly  the  semigroup  property  (1.9),  this  can  be  written  as 

S(f  +  M)  «  [Sx(t  +  k,t+  ik)S2(t  +  k,  t  +  ifc)] 

X  [S2{t+  i*,t)b’i(t+  £M)|.  ’ 


When  viewed  in  this  way  it  is  apparent  that  second  order  accuracy  may  also  be  retained 
by  using  a  splitting  of  the  form  (1.19)  but  reversing  the  order  of  S%  and  S2  in  alternate 
time  steps. 

Strang[49]  proves  that  this  splitting  is  second  order  accurate  on  a  general  nonlinear 
problem.  This  proof  is  repealed  in  Section  2.3,  where  we  also  compute  a  general 
expression  for  the  error  in  the  splitting. 

Once  the  appropriate  splitting  of  the  exact  solution  operator  has  been  chosen,  the 
time-split  method  results  from  replacing  the  exact  solution  operators  Si(Jfc)  and  S2(Jfe) 
by  approximations  Qi(k)  and  Qa(k).  A  numerical  method  based  on  the  splitting  (1.21) 
would  thus  be 

Uy-1  =  Ql{k/2)Qi(k)Qi(k/2)U^  .  (1.23) 

In  practice  Un+1  is  computed  via  the  sequence 

U'm  =  Qi(*/2)t/ 1 

UZ  =  <M*)*C  (1.24) 

uz+l  =  Qi(*/2)£C 

where  we  have  introduced  nonphysical  intermediate  solutions  U *  and  U** .  When  several 
steps  of  (1.23)  arc  applied  successively  the  adjacent  Q\[k/2)  operators  can  be  combined 
into  Qi(k),  and  the  half-step  operators  need  only  be  applied  at  the  beginning  and 
immediately  before  printout,  i.e., 

=  Qi(*/2)g2(*)Qi(fc)--Qi(fc)e2(fc)Qi(fc/2)t/2.. 

When  the  original  problem  is  split  into  more  than  two  pieces,  say 

/l(tt)  —  7fi(u)  !  ^a(w)  +  "  •  +>fp(w)» 

the  following  splitting  is  second  order  accurate: 

S(k)  »  St(*/2)S2(*/2).  •  •Sp-i(k/2)Sp(k)Sp-i(k/2y  ■  •Sa(fc/2)S1(*/2). 

This  is  easily  proved  by  induction  (see  Cottlicb[23]). 

1.4.  Hyperbolic  splittings  with  disparate  wave  speeds. 

This  thesis  is  mainly  an  investigation  into  the  applicability  of  time-split  methods  to 
pure  hyperbolic  systems  whose  solutions  consist  of  waves  traveling  at  disparate  speeds. 
Consider  the  one-dimensional  hyperbolic  constant  coefficient  system 


ut  =  Au,.  (1.25) 

The  r  X  r  matrix  A  is  assumed  to  be  diagonalizablc  with  real  eigenvalues  Mi>  l*2>  •  •  •,  Mr- 
If  X  is  the  matrix  of  right  eigenvectors  of  A,  then 

A  =  XMX~l 
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where  M  =  diag(j*i,/ia ,...,Pr)  is  a  diagonal  matrix.  The  solution  to  (1.25)  with  initial 
conditions 

«(*.0)  =  /(*) 


is  given  by 


u(z,  t)  =  exp(£/lds)ti(z,  0) 

=A-exp (tMdjX-'ftx). 


Set  v(x)  =  X~lf{x).  Then 


u(x,t)  =  X 


vi(x  +  tfil] 
vt(x  +  tp a) 

yn(x  +  tpr)  J 


In  general  each  component  tiy(z,t)  of  the  solution  is  a  linear  combination  of  waves 
traveling  at  the  various  speeds  Pi,P3,...,Pr‘  Eigenvalues  pj  with  large  amplitude  give 
rise  to  Fast  waves,  those  with  small  amplitude,  to  slow  waves. 

Suppose  now  that  the  eigenvalues  are  ordered  by  magnitude,  and  that  some  of  them 
are  much  larger  than  others: 


ImiI  <  |/*a|  <  •  •  •  <  <  K+il  <  •  •  •  <  |Mr|-  (1-26) 


Now  consider  the  use  of  a  finite  difference  scheme  for  solving  (1.25).  Throughout  this 
thesis  we  will  restrict  our  attention  to  the  Lax-Wendroff  method  for  hyperbolic  problems, 
both  in  computational  examples  and  in  some  of  the  theory  (for  example  in  Section  2.5). 
The  same  sort  of  analysis  can  be  applied  to  other  methods  with  similar  results,  but  it 
sccins  most  instructive  to  concentrate  on  one  particular  method. 

The  Lax-WendrolT  method,  like  all  explicit  methods,  is  only  conditionally  stable. 
This  places  a  restriction  on  the  size  of  the  time  step  that  can  be  used.  For  Lax-Wendroff 
this  stability  condition  is 

±p(A)  <  1.  (1.27) 

where  p(A)  =  |/ir|  is  the  spectral  radius  of  A.  The  fastest  waves  thus  dictate  the  size  of 
the  timestep  that  can  be  taken.  Accuracy  considerations  also  influence  the  size  of  the 
timestep.  In  fact  the  fastest  waves  arc  computed  most  efficiently  (in  the  sense  that  the 
least  work  is  required  to  achieve  a  given  accuracy)  if  the  mesh  ratio  k/h  rw  1  / p(A)  is 
used.  This  will  be  shown  in  Section  2.5. 

Slow  waves,  on  the  other  hand,  can  be  accurately  (and  more  efficiently)  computed 
vising  much  larger  timesteps.  The  question  is  whether  a  split  method  can  be  used  to 
compute  accurate  overall  solutions  more  efficiently. 

If  the  matrix  A  is  diagonal,  then  the  system  decouples  into  r  separate  scalar  equa¬ 
tions,  each  of  which  can  be  solved  independently  using  the  appropriate  mesh  ratio.  More 
generally,  we  can  split  the  matrix  A  into  pieces  A,  and  A /  corresponding  to  slow  waves 
and  fast  waves, 

A,  =  X  MtX~l,  Af  =  XMfX~l  (1.28) 
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where 


M,  =  diag(/*i  0, ...,  0) 

Mf  —  diag(0, . . 0,  /Vm  , . . Mr  )• 


This  essentially  decouples  the  system  into  slow  and  fast  parts.  Since  the  matrices  A ,  and 
Af  commute,  splittings  or  the  form  (1.10)  or  (1.21)  are  exact  and  nearly  optimal  mesh 
ratios  can  be  used  for  each  part. 

Realistic  problems  can  never  be  split  so  easily.  For  variable  coefficient  or  quasilinear 
systems  there  will  almost  always  be  a  splitting  error  to  contend  with.  It  is  also  generally 
undesirable  or  even  impossible  to  use  a  splitting  of  the  form  (1.28),  since  the  eigenvectors 
are  themselves  variable. 

However,  it  is  not  necessary  to  split  by  characteristic  variables  as  in  (1.28),  and 
the  time-split  method  is  often  advantageous  even  when  the  splitting  error  is  nonsero. 
Suppose,  for  example,  that  r  is  large  but  that  the  matrix  A  has  only  a  few  large 
eigenvalues.  It  may  be  the  case  that  relatively  few  elements  of  A  contribute  to  the  fast 
waves.  We  could  then  split  A  as  A  =  Aj  +  A,  in  such  a  way  that  Af  is  sparse  compared 
to  A  while  A ,  has  only  small  eigenvalues.  Because  of  its  sparsity,  taking  small  timesteps 
on  Af  requires  less  work  than  taking  small  timesteps  with  the  full  matrix  A.  The  matrix 
A,  can  be  handled  more  efficiently  using  larger  timesteps.  We  could  thus  consider  u~ing 
a  scheme  of  the  form 

tC+1  =  Qf{k/2)Q.(k)Qf(k/2)U^  (1.M) 


with 


Q.(k)  =  LW(A.,k) 

Q,(k/2)  =  (LW(A,,k/m)r /a 


(1.30) 


for  some  even  integer  m.  The  accuracy  and  the  efficiency  of  such  a  method  relative  to 
an  unsplit  method,  say  LW(A,  k/m),  depends  greatly  on  the  nature  of  the  splitting  error. 
This  will  be  studied  in  detail  in  Chapter  2. 


Example  1.1.  An  interesting  model  system  for  problems  of  this  form  is  a  block 
triangular  system  with 

\An  ^ia| 

0  A22J’ 


A  = 


(1.31) 


Suppose  the  eigenvalues  of  An  are  large  relative  to  those  of  Aaa  and  consider  the  splitting 


Ar  — 


An  0 

A  —  f« 

1  A12 

.  0  0. 

1  A22 

(1,32) 


The  effectiveness  of  the  split  method  depends  greatly  on  the  coupling  A 12  between  the 
different  time  scales.  This  is  analyzed  in  Section  2.7.  In  Section  2.8  we  present  a  simple 
procedure  for  changing  variables  to  reduce  the  coupling. 

Perturbed  problems.  The  time-split  method  can  also  be  very  effective  on  equations 
which  arc  small  perturbations  of  some  equation  for  which  the  exact  solution  operator  is 
known.  We  will  refer  to  such  problems  as  pertrubed  problems.  For  example,  consider  a 
variable  coefficient  problem  in  which  the  coefficients  have  large  mean  values  and  small 
variation.  It  may  then  be  possible  to  split  off  a  constant  coefficient  problem  u(  =  A/t», 
that  can  be  solved  exactly,  leaving  behind  the  small  perturbations  for  A,.  We  can  then 
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use  (1.29)  and  take  Q/(k/ 2)  =  cxp(^kAfdx)  with  no  error.  Since  the  dominant  part 
of  the  operator  is  being  handled  exactly,  substantial  increases  in  efficiency  are  possible. 
The  one-dimensional  shallow  water  equations  introduced  in  the  next  section  are  of  this 
form. 

More  generally  we  may  divide  the  computational  domain  into  subintervals  and  split 
out  a  different  constant  matrix  A/  on  each  subinterval.  This  might  be  appropriate  if  the 
coefficients  are  slowly  (but  widely)  varying  so  that  perturbations  about  the  local  mean 
value  are  small.  In  this  case  A/  would  be  piecewise  constant.  Alternatively  we  can  view 
this  as  a  hybrid  method  in  which  a  different  scheme  m  used  on  each  subinterval. 

In  other  cases  the  matrix  A /  may  be  variable,  but  of  a  special  form  such  that  the 
problem  ut  —  A/t»,  can  be  solved  exactly. 

Wc  continue  to  use  the  “fast”  and  “slow”  notation  even  though  for  such  perturbed 
problems  all  of  the  eigenvalues  of  A  may  be  roughly  the  same  site.  Nevertheless  in  the 
splitting  A  —  Af  +  A,  we  assume  that  A/  has  eigenvalues  much  larger  than  those  of  Aa, 
and  so  the  subproblem  ut  =  A/u,  has  waves  which  arc  fast  relative  to  those  occurring 
in  the  subproblcm  u t  =  A, t»«. 

Example  1.2.  A  simple  example  is  the  scalar  problem 

ut  =  (1  +  a(x))us  (1.33) 

where  |a(z)|  <  1  with  the  splitting 


A/  =  1,  A,  —  a(x). 

Take  k  =  2 ph  for  some  integer  p.  The  operator  exp{ \kAfdx)  is  known  exactly: 

exp(%kAfdx)u[x,t)  =  exp(phdx)u(x,  t)  ~  u(x  +  ph,  t). 

If  Lax-Wendroff  is  used  for  the  remaining  subproblem  ut  =  a(z)ux  then  the  method 
(1.24)  becomes 


LW{ct(x),  k)U*m 

Urn  +  P«(zm)(CC+ 1  -  fC— i)  +  P8«(*m)  {(a(xm+1)  +  a(xm))(tC+l  -  u'm) 

-  (a(x«.)  -  a{xm-i))[U*m  -  £/„_,)} 

=  uZ+p. 

Notice  that  even  though  this  is  a  scalar  problem,  the  operators  dx  and  a(x)dx  do  not 
commute  and  so  the  Strang  splitting  must  be  used  to  achieve  a  second  order  method. 
This  sequence  is  shown  schematically  in  Figure  1.3  for  p  =  3. 

Eliminating  the  intermediate  solutions  U  and  IJ  ,  we  can  rewrite  this  as  a  one-step 
method: 

tC+I  =  f'm+2p  +  P“(*m+p)(tf  m+2p+l  ~  ^m+Sp-t)  +  P*«(*m+ p) 

X  [(«(xm+p+l)  +  a(*m+p))(f/m+2p+l  —  ^m+8p)  (l'*M) 

—  («(*m+p)  +  a(*m+p— j))(^m+2p  —  ^m-t-Sp-l)]1 
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FIG.  1.3.  Schematic  diagram  of  the  method  (1.34)  in  split  form. 


FlG.  1.4  Stencil  for  the  method  (1.34)  viewed  as  a  one-step  method  which 
approximately  follows  the  characteristic  of  the  problem  (1.33)  (shown,  e.g., 
by  the  dotted  line).  Note  that  values  of  <*(*)  are  used  from  near  the  middle 
of  the  interval. 


The  stencil  for  this  method  is  shown  in  Figure  1.4.  The  value  1/JJ+1  is  determined  by 
the  values  of  Um  at  zm+*,_i,  »nd  xm+ap+i.  This  scheme  can  be  interpreted  as 

a  “skewed  Lax-WendrofT"  method  whose  stencil  approximately  follows  the  characteristic 
of  the  equation,  which  has  slope  —(1  +  <*(*))  at  each  point  a.  The  value  of  u(zm,t«+i) 
should  thus  be  equal  to  the  value  of  for  some  point  i  near  xm+tp.  The  exact 

location  depends  on  the  values  of  a(z)  for  all  z  between  zm  and  z.  We  thus  expect  such  a 
skewed  method*  to  be  quite  good  if  a(z)  is  small.  Just  how  good  it  is  depends  on  the  else 
of  a(z)  and  also  on  how  rapidly  a(z)  varies.  Note  that  in  the  split  method  (1.34)  only 
values  of  o(z)  near  the  middle  of  this  interval  are  used.  It  turns  out  that  the  splitting 
error  for  this  problem  depends  on  derivatives  of  a(z).  As  we  will  sec  in  Chapter  2,  when 
a(z)  is  rapidly  varying  it  is  most  efficient  to  use  small  values  of  p,  but  the  resulting 
method  is  still  more  efficient  than  using  Lax-Wendroff  on  the  unsplit  problem. 


1.5.  The  shallow  water  equations. 

Throughout  this  thesis  the  shallow  water  equations  will  be  used  as  an  example  to 
illustrate  the  general  theory  being  developed.  The  theory  applies  to  this  system  in  a  fairly 
straightforward  but  nontrivial  manner,  and  thus  studying  these  equations  provides  some 
insight  into  the  issues  which  arise  when  splitting  methods  are  applied  to  other  practical 
problems. 

In  one  space  dimension  the  shallow  water  equations  are 


til _ [it  ff|[it 

h  t~~  h  u  h 


(1.35) 


These  equations  model  flow  in  a  channel  where  g  is  the  gravitational  constant,  h  is  the 
height  of  the  water  and  it  its  velocity.  This  system  can  be  symmetrised  by  setting  4>  = 
2  \/gh  to  obtain 

r _ i  r  .  i.  11..1 

(1.36) 


We  will  make  the  realistic  asumption  that  u  is  small  compared  to  tj>  and  that  variations 
in  4>  arc  small  compared  to  some  mean  value 


\4> "  I  <  o, 

|u|  <  t^o 


(137) 


with  (<1.  Moreover  wc  assume  that  tts,  <f>x  and  higher  derivatives  arc  also  O(e<f>0).  We 
split  the  system  (1.36)  by  taking 


Af  *  "[ 


®  * o/2 

A.  =  - 

*o/2  0  J’ 

(*-*o)/2 


(4>  -  ^o)/2 

it  J’ 


(1.38) 


The  eigenvalues  of  Af  arc  ±^o/2.  The  exact  solution  operator  exp(£fc/t/ds)  can  be  used 
on  the  grid  provided  that 

£-35  “M) 

for  some  integer  p  >  1.  The  matrix  A,  has  eigenvalues  tt  ±  (^  —  ^o)/2  which  are  smaller 
by  a  factor  of  e  than  those  of  Af. 


We  will  see  in  Section  2.9  that  using  a  time-split  method  on  this  problem  reduces 
the  errors  by  a  factor  of  t  {using  the  same  amount  of  work).  The  method  is  stable  for 
the  froien-coefficient  problem,  as  seen  in  Section  3.5,  and  in  practice  is  stable  for  the 
nonlinear  system.  The  proper  specifications  of  boundary  conditions  for  this  problem  is 
discussed  in  Section  4.5. 

1.8.  A  brief  history  of  splitting  methods. 

Now  that  the  basic  form  of  the  time-split  method  and  a  wide  variety  of  possible 
splittings  have  been  introduced,  we  pause  briefly  to  review  some  of  the  extensive  work 
that  has  been  done  on  splitting  methods.  This  survey  is  far  from  complete,  but  it  provides 
some  historical  perspective  and  references,  particularly  to  the  sources  which  have  had 
the  most  impact  on  this  thesis. 

Splitting  methods  have  been  most  extensively  studied  in  the  context  of  spatial 
splittings  of  multidimensional  problems.  The  first  splittings  were  of  implicit  methods 
for  solving  parabolic  problems  and  were  also  used  as  iteration  procedures  for  solving 
steady-state  elliptic  problems. 

The  locally  onc-dimensional  methods  were  developed  primarily  by  D*Yakonov[ll], 
Marchuk[38],  Samarskii[48]  and  Yanenko(55].  The  basic  form  of  such  methods  has 
already  been  indicated  in  Section  1.2.  For  the  heat  equation  (1.1)  using  Crank- Nicolson, 
for  example,  the  scheme  is 


(1  -  **D+,Z>_,)tC  =  (1  +  i *D+,D_,)t/«  .  40) 

(1  -  }kD+yD-.v)l/Z+i  =  (1  + 

This  clearly  fits  into  the  genera!  framework  introduced  in  the  Section  1.3  with  the  splitting 

A,(u)  =  u„,  **{«)  =  *„.  (1.41) 


The  LOD  method,  however,  was  not  the  first  such  splitting  method  to  be  used. 
In  the  mid  1950’s  the  Alternating  Direction  Implicit  (ADI)  method  was  introduced  by 
Pcaccman  &  Racliford[45]  and  Douglas[6].  On  the  equation  (1.1)  this  method,  known  as 
the  Peaccman-Rachford  method,  has  the  form 

(1  -  l,kD+,D.,)U‘m  =  (1  +  \kD+,D-,)UX 
(1  -  IkD+yD-jUX'  -  (1  +  i*D+./?_.)lC 

The  philosophy  behind  the  ADI  method  is  somewhat  different  from  that  behind  the  LOD 
method.  Each  equation  of  (1.42)  is,  by  itself,  a  first  order  accurate  scheme  for  solving 
the  original  equation  (1.1)  on  a  timestep  of  length  k/2.  The  combination  provides  a 
second  order  accurate  solution  on  a  step  of  length  k.  In  some  sense  it  is  thus  a  more 
natural  approach  to  solving  the  problem  than  the  LOD  method,  since  the  individual 
equations  composing  (1.40)  do  not,  by  themselves,  provide  consistent  approximations  to 
the  original  system.  On  the  other  hand,  the  LOD  method  can  be  viewed  more  naturally 
as  a  time-split  method  of  the  form  discussed  in  Section  1.3,  since  each  equation  of  (1.40) 
is  a  second  order  accurate  approximation  to  one  of  the  subproblcms  determined  by  (1.41) 
on  a  timestep  of  length  k. 
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In  Tact,  the  Peaccman-Rachford  method  can  also  be  viewed  as  a  time-split  method 
of  this  form,  but  with  a  different  splitting.  Instead  of  (1.41)  suppose  we  split  the  operator 
A(v)  —  u,,  +  Uyy  as 

^i(u)  —  i(ur*  +  ii»v)  + 

^lM  =  K«,.  +  “  uvnv)‘ 

Then  the  equations  of  (1.42)  are  second  order  accurate  approximations  to  u»  =  4i(«) 
and  tt*  —  ^j(u)  on  Umesteps  of  length  k. 

There  are  many  other  ways  of  relating  the  ADI  and  LOD  methods,  see  for  example 
Gourlay  &  Mitchell(25]  or  Morris  &  Gourlay(42].  One  advantage  of  viewing  ADI  aa 
a  time-split  method  is  that,  in  some  cases,  appropriate  boundary  conditions  for  the 
intermediate  solution  U *  can  then  be  easily  determined  using  the  general  procedure 
described  in  Chapter  4.  This  is  discussed  in  Section  5.4. 

Numerous  variations  on  the  Peaceman-Rachford  method  have  been  proposed  over 
the  years,  for  example  by  Douglas  &  Rachford[8],  Fairweather  &  Mitchell[19],  Douglas 
&  Gunn(7],  and  D'Yakonov(l2|.  The  last  of  these  is  particularly  interesting  since  it 
is  based  on  approximate  factorisation,  an  approach  that  is  currently  quite  popular 
in  computational  fluid  dynamics.  D'Yakanov's  method,  which  he  calls  the  method  of 
disintegrating  operators,  results  from  the  approximations 

1  —  £fc(D+*D_x  +  D+yD_y)  »  (1  -  |fcD+xD_x)(l  —  JfcD+vD_y), 

1  +  JkjD+jD-*  +  D+yD-y)  w  (1  +  ^fcD+xD_z)(l  +  jfcD+yD_y). 

Each  of  these  has  an  error  $JtaD+xD_xD+yD_y.  When  both  approximations  arc  used 
in  (1.16),  the  rcsuling  error  is  0(*3).  This  leads  to  the  split  method 

(1  -  lkD+sD.s)U'm  =  (1  +  **D+XZ>_X)(  1  +  ^D+yD_y)l/- 
(1  -  i*D+yD-y)f/-+‘  =Ul 

which  can  also  be  viewed  ns  a  time-split  method  with  the  splitting 

^|{«)  =  U„  +  }Uyy  -  £  kllyyyy 

Aa(ll)  =  Jtiyy  +  ^kUyyyy. 

A  great  deal  of  work  has  gone  into  the  proper  specification  of  intermediate  bound¬ 
ary  conditions  for  such  splitting  methods.  Sec,  for  example,  Lawson  &  Morris[84], 
Fairweather  &  Milchcll(19],  or  D’Yakonov[13].  General  discussions  of  splitting  methods 
for  parabolic  problems  can  be  found  in  many  places,  including  Yanenko|56),  Marchuk[40], 
and  Mitchell  &  Griinihs(41], 

As  opposed  to  parabolic  problems,  many  hyperbolic  systems  of  equations  arc  solved 
using  explicit  methods.  As  we  saw  for  the  one-dimensional  system  (1.25),  the  stability 
limit  frequently  allows  timesteps  that  are  reasonable  from  the  standpoint  of  efficiency, 
and  so  there  is  no  need  to  use  implicit  methods.  In  more  space  dimensions,  however, 
the  stability  limit  is  often  severely  reduced.  (For  example,  the  stability  limit  for  2D  Lax- 
Wendroff  on  (1.3)  is  Xma x(p(A),p(n))  <  l/>/8.)  Strang[49j  showed  that  if  the  locally 
one-dimensional  method  is  used  on  (1.3),  then  the  stability  limit  is  more  reasonable  (for 


Lax-Wondroff,  X  max(p(A),  p(B))  <  1).  Thus  the  LOD  method  has  the  use,  for  explicit 
methods,  of  increasing  the  stability  limit. 

Implicit  methods  are  often  used  for  certain  classes  of  hyperbolic  systems.  Recall  that 
the  timestep  for  an  explicit  method  is  restricted  by  the  fastest  wave  speed.  For  certain 
systems  of  equations  with  disparate  wave  speeds  the  physically  meaningful  solutions  con¬ 
tain  no  fast-wave  components,  or  at  least  the  fast  waves  have  small  amplitude  compared 
to  the  slower  waves.  For  an  explicit  method  applied  to  such  a  problem,  stability  con¬ 
siderations  limit  the  timestep  to  a  value  much  smaller  than  that  required  for  accuracy. 
For  this  reason,  implicit  methods  are  frequently  used  instead.  In  more  than  one  Bpace 
dimension  LOD,  ADI  or  approximate  factorization  methods  again  prove  useful. 

Such  problems  arise,  for  example,  in  modeling  atmospheric  flows.  The  simplest 
such  system  is  the  two-dimensional  shallow  water  equations.  The  general  solution  to 
these  equations  includes  both  fast  “gravity  waves”  and  much  slower  “Rossby  waves”.  In 
practice,  however,  the  gravity  waves  contain  little  energy  and,  it  is  thought,  have  little 
effect  on  the  weather.  Gustafsson[29]  has  studied  an  ADI  method  for  this  problem. 

Another  approach  to  such  problems  has  been  taken  by  Kreiss[32],[33]  and  Browning, 
Kasahara  &  Krciss(3].  They  properly  prepare  the  data  so  that  fast  wave  components  are 
eliminated.  Majda[37]  has  considered  using  filters  to  suppress  the  fast  waves  in  the  same 
context. 

Approximate  factorization  methods  for  hyperbolic  problems  have  been  studied  by 
Warming  &  Beam[54],  primarily  for  the  Euler  equations  of  gas  dynamics  and  for  mixed 
hyperbolic- parabolic  problems  such  as  the  Navicr-Stokes  equations.  Again  they  are 
dealing  with  problems  where  the  fast  waves  have  little  effect  on  the  solutions  of  interest. 

Another  possible  approach  for  such  problems  is  to  split  the  coefficients  into  fast  and 
slow  terms  and  use  an  implicit  method  only  on  the  fast  part.  This  can  be  quite  efficient  if, 
for  example,  the  fast  part  is  sparse.  The  splitting  between  implicit  and  explicit  methods 
can  be  effected  in  various  ways.  For  the  problem  «*  =  Auz  =  (A/  +  A,)ttz,  a  time- 
split  method  of  the  form  (1.20)  could  be  used  with  Qf{k/2)  implicit  and  Qa{k)  explicit. 
Alternatively,  one-step  methods  can  be  derived  that  arc  implicit  only  in  A/.  For  example, 
the  trapezoidal  formula  and  leapfrog  can  be  combined  into  the  hybrid  method 

(7  -  kAfD0)U •»+*  =  (7  +  kAfDoWZT1  +  2kA,D0U -  .  (1.44) 

Such  methods  arc  called  semi-implicit  methods  or  explicit- implicit  methods.  Elvius  & 
Sundstrom[16]  have  analyzed  the  two-dimensional  analogue  of  (1.44)  for  the  shallow 
water  equations.  Harlow  &  Amsden(31]  have  applied  a  similar  method  to  the  Euler 
equations. 

The  idea  of  using  different  timesteps  on  various  parts  of  the  system  has  been  used  in 
one  form  or  another  by  several  authors,  including  Engquist,  Gustafsson  &  Vreeburg(17|, 
Gadd[20],  and  Turkel  Sr.  %was[52]. 

Many  nonlinear  hyperbolic  systems  have  solutions  involving  shock  waves  discontinu¬ 
ous  solutions  which  can  arise  even  from  smooth  initial  data.  For  such  problems  a  wide 
variety  of  special  methods  have  been  devised.  Often  these  methods  are  directly  applicable 
only  in  one  space  dimension.  For  higher  dimensional  problems,  LOD  splittings  arc  agaii. 
frequently  used.  Since  the  solutions  arc  not  smooth,  splittings  are  more  difficult  to 
analyze  in  this  context.  Crandall  &  Majda[5]  have  proved  that  the  splittings  (1.19)  and 
(1.21)  do  give  convergent  methods  when  applied  to  scalar  conservation  laws. 
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For  mixed  problems  such  as  the  Navier-Stokes  equations  the  time-split  method  has 
been  used  to  split  between  hyperbolic  and  parabolic  parts.  Abarbanel  &  Gottlicb(l)  split 
the  full  three-dimensional  Navier-Stokes  equations  into  nine  pieces — the  hyperbolic  and 
parabolic  terms  in  each  space  dimension  and  three  cross-derivative  terms.  They  then 
use  the  time-split  method  to  derive  an  explicit  method  with  good  stability  properties. 
MacCormack[35j[36],  Strikwerda[50],  and  Dwoycr  &  Thames[10]  have  studied  similar 
methods. 

Approximate  factorisation  methods  for  this  problem  have  been  proposed  by  Beam 
and  Warming[2](54].  This  approach  appears  to  have  certain  advantages  in  steady  state 
calculations.  The  numerical  steady  state  is  independent  of  the  timestep  k  used  to  compute 
it,  and  the  calculations  can  be  performed  in  an  “increment  form”  that  is  computationally 
efficient. 

Convection-diffusion  equations  similar  to  (1.4)  arise  in  transport  problems  that  in¬ 
clude  diffusion,  for  example  in  multi-phase  miscible  flow  or  in  modeling  heat  flow  in  a 
moving  material.  When  the  problem  is  convection  dominated  (u  <C  |c|  in  (1.4)),  the 
propagation  of  sharp  fronts  is  often  of  interest.  These  are  difficult  to  handle  numerically. 
It  is  often  advantageous  to  again  split  between  the  hyperbolic  and  parabolic  parts  and 
handle  the  hyperbolic  part  using  characteristics.  This  is  studied  in  Section  5.3.  Similar 
methods  have  been  proposed  by  MacCormack[36]  and  Douglas  &  Russell[9).  Another 
possibility  is  to  use  the  Unite  element  method  for  the  parabolic  parl[9][15]{47][53]. 

1.7.  Outline  and  summary  of  results. 

There  are  three  main  issues  to  be  dealt  with  when  considering  the  use  of  a  time-split 
method  for  any  differential  equation.  These  may  be  summarized  as  efficiency,  stability, 
and  the  proper  choice  of  boundary  conditions.  These  are,  of  course,  major  issues  in  the 
choice  of  any  Finite  difference  scheme,  but  the  use  of  time-split  methods  introduces  new 
complications  into  each  area. 

Efficiency.  The  first  quantity  to  compute  in  the  analysis  of  any  finite  difference 
scheme  is  its  truncation  error.  In  Section  2.2  we  show  that  Tor  the  time-split  method  the 
truncation  error  is  simply  the  sum  of  the  splitting  error  and  the  truncation  errors  for 
the  approximate  solution  operators  Qi  and  Q?  (plus  higher  order  terms).  The  splitting 
error  thus  plays  a  fundamental  role  and  techniques  for  computing  this  error  for  general 
splittings  arc  discussed  in  Section  2.3. 

In  comparing  methods,  however,  it  is  not  in  general  sufficient  to  compare  their 
truncation  errors,  since  one  scheme  may  require  much  more  computation  than  another. 
This  is  particularly  true  where  time-split  methods  arc  concerned.  Instead  we  must 
compare  some  measure  of  the  efficiency  of  the  methods  such  as  the  amount  of  work 
required  to  achieve  a  given  accuracy. 

Since  split  methods  generally  involve  the  conjunction  of  several  different  numerical 
methods,  there  may  be  several  free  parameters,  such  as  stepsizes,  to  be  chosen.  For  the 
method  (1.30),  for  example,  we  must  choose  both  k/h  and  m.  We  can  essentially  choose 
the  mesh  ratios  for  the  two  time  scales  independently.  As  we  will  see  in  Section  2.5, 
the  optimal  choice  depends  on  the  size  of  the  splitting  error  and  is  not  always  obvious  a 
priori,  in  particular,  the  optimal  mesh  ratio  k/h  is  often  Tar  from  the  stability  limit  of 
the  method  used  on  the  slow  problem. 


Stability.  For  some  time-split  methods  applied  to  certain  problems,  the  operator 
Qi(k)Q*(k)  is  stable  whenever  the  operators  Qi(k)  and  Qt{k)  are  each  stable  operators  on 
their  own.  Unfortunately,  this  is  not  true  in  general;  the  product  of  two  stable  operators 
may  be  unstable.  An  example  of  this  is  given  in  Chapter  3. 

Of  course  the  stability  of  Qtj[k)Qa(k)  can  always  be  determined  directly  by  eliminat¬ 
ing  all  intermediate  variables  and  viewing  the  split  method  as  a  one-step  method.  How¬ 
ever,  the  resulting  method  is  generally  quite  complicated,  making  direct  analysis  difficult. 
For  this  reason  it  is  useful  to  identify  special  classes  of  problems  for  which  the  individual 
stability  of  Qi(k)  and  Qa(k)  does  guarantee  the  stability  of  Qt(k)Qa(k).  Several  such 
classes  of  hyperbolic  splittings  are  identified  in  Chapter  3. 

Boundary  conditions.  Ail  practical  calculations  are  performed  on  finite  domains. 
If  periodic  boundary  conditions  are  used  (e.g.,  u(0,  t)  =  u(l,t)  Vf  on  the  strip  0  <  x  < 
1),  then  the  same  finite  difference  scheme  can  be  used  at  all  points,  simply  by  wrapping 
around  at  the  boundaries.  Otherwise,  one  or  more  points  at  each  boundary  will  have 
to  be  determined  in  some  alternative  manner  (unless  a  one-sided  scheme  is  used).  Some 
boundary  values  will  be  provided  as  part  of  the  problem,  but  frequently  finite  difference 
approximations  require  more  boundary  conditions  than  the  original  differential  equation. 
The  remaining  boundary  values  must  be  determined  by  some  other  procedure.  A  variety 
of  techniques  arc  used  for  this  purpose,  depending  on  the  context.  The  easiest  approach  U 
often  to  extrapolate  the  interior  solution  at  time  tn+ j  oul  to  the  boundary.  Alternatively, 
one-sided  (or  lopsided)  finite  difference  schemes  can  be  used  to  compute  the  solution  at 
points  on  (or  near)  the  boundary.  At  some  boundaries  other  desirable  properties  of  the 
solutions,  such  as  nonreflection  of  outgoing  waves,  may  be  used  to  determine  the  proper 
boundary  values. 

For  time-split  methods  the  choice  of  boundary  values  is  complicated  by  the  need 
to  supply  boundary  data  for  the  intermediate  solutions,  such  as  U*.  These  solutions 
arc  obtained  not  by  solving  the  original  differential  equation  but  rather  by  solving  one 
of  the  subproblcms.  Because  of  this,  appropriate  boundary  data  for  the  intermediate 
solutions  is  never  available  directly.  [Extrapolation  from  the  interior  can  still  be  used, 
but  is  generally  undesirable  both  for  reasons  of  stability  and  accuracy. 

The  generation  of  boundary  data  for  the  intermediate  solutions  is  discussed  in 
Chapter  4.  We  describe  a  general  procedure  for  transforming  given  boundary  data  for  the 
original  equation  into  appropriate  data  for  the  intermediate  solutions.  This  procedure 
is  based  on  the  following  idea.  We  introduce  a  new  function  u  which  satisfies  the 
subproblein  that  is  actually  being  solved  in  the  first  step  of  the  splitting.  We  then  expand 
the  desired  boundary  value  for  u*  in  a  Taylor  series  about  the  initial  time  tn  at  which 
u  =  u.  Using  the  differential  equations  for  u*  and  u  we  then  reexpress  this  as  a  scries 
expansion  involving  only  the  function  u  and  its  time  derivatives  along  the  boundary. 
This  can  then  be  evaluated  from  the  given  boundary  conditions  for  «. 

fiach  of  the  next  three  chapters  is  devoted  to  one  of  these  issues.  The  emphasis 
is  on  splittings  of  hyperbolic  problems  iuto  subproblcms  with  disparate  wave  speeds,  as 
discussed  in  Section  1.3.  However,  many  of  the  techniques  used  arc  also  applicable  to 
other  splittings  of  the  form  ut  =  ^i(u)  +  ^a(u).  Whenever  possible,  the  discussion  is  in 
terms  of  the  more  general  splitting  to  facilitate  application  to  other  problems.  Hyperbolic 
splittings  are  always  used  as  concrete  examples  in  these  chapters,  and  most  of  the  specific 
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results  arc  for  such  problems.  In  particular,  the  one-dimensional  shallow  water  equations 
are  frequently  used  as  an  example. 

In  Chapter  5  we  discuss  several  other  applications  of  the  time-split  method  using  the 
theory  developed  jn  previous  chapters.  We  first  consider  applications  of  the  time-split 
method  to  hyperbolic  problems  in  two  space  dimensions.  The  main  intent  is  still  to  split 
between  different  wave  speeds,  but  in  conjunction  with  this  spatial  splittings  may  also 
be  used. 

Finally  we  consider  two  applications  of  the  theory  of  time-split  methods  to  non- 
hyperbolic  splittings.  In  Section  5.3  the  simple  convection-diffusion  equation  (1.4)  is  split 
and  solved  as  a  perturbed  problem  with  a  skewed  Crank-Nicolson  method  analogous  to 
the  skewed  Lax-Wcndroff  method  (1.34).  The  efficiency  of  this  method  can  be  analysed 
using  the  techniques  of  Chapter  2.  Intermediate  boundary  conditions  at  the  inflow  bound¬ 
ary  can  be  specified  using  the  procedure  of  Chapter  4.  When  the  diffusive  parameter  v 
in  (1.4)  is  small  the  equation  becomes  a  singular  perturbation  equation  with  a  boundary 
layer  at  the  outflow  boundary  that  causes  additional  difficulties^ 

In  Section  5.4  the  Peaceman-Rachford  method  (1.42)  is  viewed  as  a  time-split  method 
with  the  splitting  (1.43).  For  a  rectangular  region  the  boundary  condition  procedure  of 
Chapter  4  can  be  used  to  derive  appropriate  boundary  conditions  for  the  intermediate 
solution  U  .  These  arc  seen  to  agree  with  the  classical  boundary  conditions  of  Fairweather 
and  Mitchcll[19|. 

Chapters  2-4  are  essentially  independent  of  one  another  and  can  be  read  in  any 
order.  The  sections  of  Chapter  5,  which  deal  with  other  applications,  are  disjoint  from 
one  another,  but  build  upon  the  results  of  the  previous  chapters,  particularly  Chapters 
2  and  4. 


2.  Accuracy  and  efficiency 


2.1.  Introduction. 

This  chapter  begins  with  a  computation  of  the  truncation  error  for  a  general  time- 
split  method.  Neglecting  higher  order  terms,  this  is  simply  the  sum  of  the  error  com¬ 
mitted  in  splitting  the  exact  solution  operator  (the  splitting  error)  and  the  truncation 
errors  of  the  schemes  used  for  the  subproblems. 

In  Section  2.3  we  present  general  expressions  for  the  splitting  error  in  both  the  first 
order  splitting  (1.19)  and  the  Strang  splitting  (1.21).  The  splitting  error  is  explicitly 
computed  for  some  model  problems,  including  the  one-dimensional  shallow  water  equa¬ 
tions. 

For  the  type  of  splitting  with  which  this  thesis  in  most  concerned,  namely  where 
||/?2(«)||  <  tll^i(«)||  with  £  4C  1,  the  error  in  the  Strang  splitting  is  seen  to  be  0(ck3).  A 
simple  modification  of  this  splitting  is  proposed  with  0(c3k3  +  tk*)  splitting  error. 

Once  wc  arc  able  to  compute  the  splitting  error  for  specific  problems,  we  can  analyze  ' 
the  efficiency  of  the  split  method  relative  to  unsplit  methods.  It  turns  out  that  the  size  of 
the  splitting  error  greatly  affects  what  size  timesteps  should  be  used  in  the  split  method 
and  what  increase  in  efficiency  can  then  be  expe.cted.  This  analysis  is  presented  in  Section 
2.5  and  continues  in  Section  2.6  where  phase  errors  arc  computed. 

In  Section  2.7  these  results  are  interpreted  for  a  block  triangular  system  of  the  form 
considered  in  Kxainplc  1.1.  For  this  problem  (and  also  for  more  general  partitioned 
systems)  the  splitting  error  can  be  reduced  by  the  use  of  a  simple  change  of  variables. 
This  is  discussed  in  Section  2.8. 

In  Section  2.9  the  onc-dimcnsional  shallow  water  equations  arc  studied.  The  theory 
developed  in  Section  2.5  is  confirmed  numerically  Tor  this  system. 

2.2.  TVuncation  error  of  the  time-split  method. 

In  order  to  compute  the  truncation  error  for  the  time-split  method  we  first  introdupe 
the  truncation  error  operators  fi,{k)  for  the  approximate  solution  operators  Qj(k).  These 
arc  defined  by 

/-;,(*)  =  Qj(k)  -  Sj(k),  k  —  1,2. 

We  will  assume  throughout  that  Q\  and  are  at  least  second  order  accurate.  Then 
Ej(k)u  is  ()(k3)  for  smooth  u.  For  shorthand  wc  sometimes  write  Iij(k)  =  0{k3). 
Similarly,  we  introduce  the  splitting  error  operator  K,piit (k)  defined  by 


/i;„iit(*)  =  5,(*/2)5'a(ifc)S1(*/2)  -  S(k). 


This  is  also  0(k 3)  for  smooth  u. 


The  truncation  error  operator  for  the  time-split  method  is 

ETau(k)  =  Qi(k/2)Q2(k)Qi(k/2)  -  S(k). 

If  the  operators  and  A2(u)  are  linear,  this  can  be  easily  computed  in  terms  of  E\, 
E2  and  B,pih  using  the  fact  that  Sj(k)  =  I  +  O(k)  and  Qj(k)  =  I  +  O(k): 

ETa“(k)  =  (S,(fc/2)  +  Et(k/ 2))  (S2(k)  +  E2(k))  (St{k/ 2)  +  E^k/2))  -  S(k) 

=  Si(k/2)S2(k)Si(k/2)  -  S(k)  +  2E,(Jfc/2)  +  E2{k)  +  0{k4) 

=  E'Plit(k)  +  2El(k/2)  +  E2(k)  +  0(k4).  (2.1) 

If  the  operators  ^i(u)  and  ^(u)  are  nonlinear,  then  the  Q,  S,  and  E  operators  will  also 
be  nonlinear.  More  care  must  then  be  used  in  deriving  ET3u(k),  but  the  0(k3)  term  of 
the  result  is  exactly  the  same  as  above,  and  the  expression  (2.1)  holds  in  general. 

The  truncation  error  for  the  time-split  method  is  thus  seen  to  be  essentially  the  sum 
of  the  splitting  error  for  the  problem  and  the  truncation  errors  for  the  finite  difference 
operators.  This  allows  us  to  easily  compute  the  accuracy  and  investigate  the  efficiency 
of  the  time-split  method  relative  to  unsplit  schemes.  This  will  be  done  in  Section  2.5. 
First  we  must  be  able  to  compute  the  splitting  error  E»piu(A:). 


2.3.  The  splitting  error. 

We  will  first  prove  the  assertions  made  in  Chapter  1  regarding  the  accuracy  of  the 
splittings  (1.19)  and  (1.21)  when  applied  to  the  solution  operator  for  a  general  equation 
of  the  form 

«t  =  >l(«,  t).  (2.2) 


The  operator  A  may  also  depend  on  spatial  variables,  but  this  dependence  will  not 
be  explicitly  shown.  The  proofs  arc  completely  independent  of  the  number  of  space 
dimensions. 

We  denote  by  A'(u,  t)  the  total  time  derivative  of  A  assuming  u  satisfies  (2.2).  This 
is  given  by 

A'(u,  t)  =  At(n,  t)  +  Au(u,  t)u, 

=  At(u,  1)  +  Au(n,  t)A(u,  t). 


In  the  latter  form  this  depends  only  on  u  at  the  time  t  and  does  not  depend  explicitly 
on  ut.  This  is  crucial  in  the  proofs  that  follow,  where  we  will  be  switching  between 
solving  different  differential  equations,  which  means  that  time  derivatives  of  tt  become 
ambiguous. 

A  few  words  should  be  said  about  the  quantity  >!„(«,  t).  The  vector-valued  function 
A  generally  depends  both  on  u  and  on  one  or  more  spatial  derivatives  of  «.  In  one  space 
dimension,  for  example,  we  could  write  A  =  A{u,  ux,  «„,... ,f).  The  derivative  Au  is 


then  defined  aa 


„  dA  dA  _  dA  . 
Am  =  - - h  — - <3*  +  ~Z - 1 

uli  uUg  vUu 


(2.4) 


where  OA/dn,  dA/Ous,  etc.  are  ordinary  Jacobian  matrices  with  respect  to  the  ap¬ 
propriate  vectors  «,  <iz,  etc.  More  will  be  said  about  evaluating  these  expressions  later 
in  this  section. 
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Now  suppose  that  A  is  split  as  A(u,  t)  —  ^t(u,f)  +  Ai(u,  t).  One  consequence  of  (2.3) 
is  that  A'(u,t)  A\(u,t)  +  A2(u,t).  This  is  because  A'j(u,  t)  is  the  time-derivative  of  Aj 

assuming  a  satisfies  ut  =  Aj( a,  t)  rather  than  (2.2).  We  find  instead  that 

A'  =  At  +  AuA 

=  (Ae  +  An)  +  {Ai*  +  ^2«)(A  +  ^a)  ^2  5) 

=  (^it  +  AiuAi)  +  (An  +  AiuAi)  +  AiuAi  +  Ai^At 

=  A\  +  A'i  +  Ai^At  +  AiuAi- 

We  are  now  ready  to  prove  the  results  indicated  earlier,  beginning  with  a  standard 
proof  that  the  splitting  (1.19)  is  first  order  accurate  (i.e.,  that  the  local  error  is  0(k2)). 

THEOREM  2.1.  Suppose  that  u(to)  is  a  C°°  function  of  all  spatial  variables  and  that 
A,  At,  and  A%  are  smooth  functions  of  a  and  t  related  by  (1.17).  Then  the  corresponding 
solution  operators  S,  Si,  and  S2  satisfy 

Si(t0  +  k,t0)Si(t0  +  k,t0)u(to)  —  S(to  4-  k,to)u(to) 

=  5^2[>l2ii(«(^o)i  to)At(ti(io),  to)  ~  Aiu(u(to)t  to)Aii<.{to),  to))  +  0(k 3) 

as  k  — *  0. 

Proof.  We  begin  by  computing  S(to  +  k,  to)u(to).  If  a  satisfies  (2.2)  then  this  is 
simply  a(to  +  k)  and  expanding  in  a  Taylor  scries  gives 

S(to  +  A;,  to)a(to)  ==  u(to)  +  kut(to)  +  ^k2utt(to)  +  0(k 3)  . 

=  u(t0)  +  kA  +  ik2A'  +  O(k3).  [  '  ’ 

Here  and  below,  when  no  arguments  are  shown  for  A  we  mean  /l(«(to).to)  (similarly  for 
A 1  and  A2).  We  now  compute  the  solution  using  the  split  operator.  Alter  the  first  step 
we  have 

5, (t0  +  k,  fo)u(<o)  —  a(to)  +  kA  1  +  ^k2A\  +  0(k3). 

Set  u  =  5i(to+  k,t0)u(t0).  Then 

Si(to  +  k,  to)u  =  a  +  kA%(u  ,  to;  +  gA2^2(a  ,to)  +  0(A3) 

=  a  +  k[Az(u(to),  to)  +  ^2«(a(to)>  to)(a  —  a(to))  +  0(A2)) 

+  A2[>i,2(a(fo),  t0)  +  O(A)]  +  0(k3) 

=  [u(t0)  +  kAi  +  i2k2A\+O(k3)} 

+  k[A2  +  A2»(kAl  +  0(k2))  +  0(k2)} 

+  ^k2[A'2  +  0(k)\  +  0(k3) 

=  a(t0)  +  k(Ai  +  At)  +  \k2(A\  +  2A2»Al  +  A'2)  +  0(k3). 

Using  (2.5)  and  (2.7)  wc  find  that 

S2(to  +  k,  fo)5i(io  4-  k,  to)a(to)  —  S(to  +  k,  to)a(to) 

=  lk2(A2uAi-AluA2)  +  0(k3). 

This  proves  the  theorem.  | 
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Example  2.1.  The  formula  (2.6)  can  be  used  to  compute  the  0(A2)  term  of  the 
splitting  error  for  any  particular  problem.  Consider,  for  example,  the  problems 

(a)  ut  =  Aus  +  Buv  with  Aj(u)  =  Aux,  At(u)  —  Buv  and  A  and  B  constant. 
The  solution  operators  are  simply  exponentials,  so  the  splitting  error  can  be  computed 
directly  as  in  (1.20).  We  get  the  same  result  by  (2.6)  since  Ai*  =  Adx  and  At%  —  Bdv. 

(b)  ut  =  (1 4- a(z))ux  with  /fi(u)  =  ux  and  ^(u)  =  a(z)us.  From  (2.6)  the  splitting 
error  is 

^ki[a(x)dxus  -  3*(a(z)tt,)]  +  0(AS) 

=  -  £A2a'(z) «*  +  0(A3). 

For  this  problem  the  solution  operators  are  again  exponentials,  Si(k)  =  exp(Adx)  and 
52(A)  =  exp(Aa(z)5x),  so  this  can  also  be  checked  directly. 

(c)  ut  =  ( c+u)ux  with  c  constant  and  u  scalar.  Take  Ai(u)  —  cur  and  At(u)  =  uux. 
Then  A\xAt  =  AtuAt  ==  c(ux  4-  u«„)  and  the  0(k2)  term  in  the  splitting  error  is  zero. 
In  fact,  for  this  problem  the  splitting  error  is  identically  zero.  This  is  intuitively  clear 
since  solving  the  subproblem  ut  =  cux  simply  translates  the  solution  in  z.  The  remaining 
subproblem  ut  =  uux  does  not  depend  explicitly  on  z,  and  so  solving  this  problem  and 
shifting  the  result  is  equivalent  to  solving  the  original  problem. 

Note  that  if  u  is  a  vector  this  is  no  longer  true,  since  in  general  different  eigen- 
components  of  u  propagate  at  different  speeds  and  hence  move  relative  to  one  another. 
The  splitting  error  for  a  system  of  equations  ut  =  \Af  4- -/4#(u))ux  will  be  computed  later 
in  this  section. 

The  next  theorem  asserts  that  the  Strang  splitting  (1.21)  is  in  general  second  order  - 
accurate. 

THEOREM  2.2.  (Strang[ 49])  Suppose  that  u(to)  is  a  C°°  function  of  all  spatial 
variables  and  that  A,  At,  and  At  are  smooth  functions  of  u  and  t  related  by  (1.17).  Then 
the  corresponding  solution  operators  S,  Si,  and  St  satisfy 

Si(t0  4-  k,  t0  +  fk)St(to  4-  k,to)Si(to  +  gA,  t0)u(to)  —  S(to  4-  A,  to)u(to)  =  C(A3) 
as  A  — >  0. 

Proof.  Proceeding  as  in  the  proof  of  Theorem  2.1, 

S\  (to  +  jA,fo)w(to)  =  «(*o)  +  5AA1  4-  fak2A\  4-  0(A3). 

Again  denote  this  by  «*.  Then 

■^(fo  +  A,  to)u  =  «  +  k Ailit  , to)  4-  jAayJj(tt  ,  to)4-C(A3) 

—  u(to)  4-  k(£Ai  4-  At)  +  yA2(j^j  4-  At*A\  4-  Aj)  4-  0(k3). 

Call  this  quantity  u**.  Then 

•**1  (to  +  A, to  4-  j A)«  =  u  4- ^Ayli(u  ,to  4*  ^A)  4- |Aaylj(u  4- ^A)  4- 0(A3). 

Expanding  A 1  and  A\  in  both  u  and  t  about  (u(f0),  to)  and  collecting  terms,  we  find  that 
i’.(t0  4-  A,  to  4-  gA)5a(fo  +  A,to)6’i(fo  4-  jA,  to) 

=  u(to)  +  k(Ai  +  At)  4-  2 A^[ A\  4-  4-  AiuAi)  4-  A 3 

4-  Ai„At  4-  ^2u^i]  +  0(A3) 

=  u(t0)  +  kA  +  JA2/!'  4-  0(A3) 

in  view  of  (2.3)  and  (2.5).  Comparing  this  with  (2.7)  shows  that  the  error  is  indeed 

0(A3).  | 
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Keeping  the  0(k3)  term  everywhere  in  the  proof  would  have  given  us  a  formula 
analogous  to  (2.6)  for  the  k 3  term  of  the  error.  In  general  this  is  quite  complicated.  For 
the  relatively  simple  autonomous  case  where  A  is  a  function  of  u  alone,  the  splitting  error 
is  found  to  be 


— .^k3[^AiuAi)^dt  —  +  %iAa*Ai)nAt 

—  %.AtuAa)uAi  +  (AauAi)uAa  —  %AiuAa)uAa]  +  0(*4)- 


(2.8) 


Example  2.2.  The  errors  in  the  Strang  splitting  for  the  problems  considered  in 
Example  2.1  are  relatively  easy  to  compute: 

(a)  ut  =  Aux  +  Buy  with  A  and  B  constant.  By  expanding  the  exponential  solution 
operators  the  splitting  error  is  seen  to  be 

-  **3((KS  -  hABA  +  \BA3)d\dy 

-  ( %B2A  -  BAB  +  %AB*)dsdl]u(t0)  +  0(fc4).  V  ‘  ' 

The  splitting  error  is  zero  only  if  A  and  B  commute. 

(b)  ut  =  (1  +  a(x))ux.  Again  expanding  the  exponential  solution  operators  shown 
that  the  splitting  error  is 


-  ***((i+o(«))a'(x)  -  (a'(*))8Mt0)  +  0[k<). 


A  higher  order  splitting.  The  fact  that  the  Strang  splitting  is  second  order 
accurate  can  be  seen  more  directly  by  viewing  the  Strang  splitting,  as  in  (1-22),  as  two 
applications  of  the  first  order  splitting  with  Si  and  Sa  applied  in  the  opposite  order  in 
the  second  application.  By  Theorem  2.1  the  truncation  error  in  the  first  step  is 

£*2(yt2»(«(*o),  t0)AMt0),  to)  -  AMto),  t0)Aa[u(t0),  to))  +  0(k3)  (2.10) 

and  in  the  second  step: 

ffc2Mtu(u  ,to  +  j^)^2(u  i  to  +  jfc)  ~  ^2»(u  >to+5^Mi(w  i  to  +  5*))  +  0(fc3).  (2.11) 

The  full-step  truncation  error  can  be  shown  to  be  simply  the  sum  of  (2.10)  and  (2.11) 
plus  0(A:3)  terms.  Expanding  (2.11)  about  (u(fo),  to)  and  adding  (2.10),  the  0(k 2)  terms 
cancel  and  hence  the  Strang  splitting  is  0{k3)  accurate.  This  cancellation  occurs  because 
the  0(k 2)  term  of  (2.6)  is  skew-symmetric  in  the  variables  Ai  and  Aa- 

For  the  type  of  problem  we  are  considering  here,  where  ||>l2(u)||  <  f||>i|(ti)||  and 
similarly  for  their  derivatives,  a  similar  trick  can  be  applied  to  the  Strang  splitting  to 
increase  the  accuracy  even  further.  The  0(fc3)  term  of  the  splitting  error  (2.8)  is  generally 
dominated  by  the  first  three  terms,  which  contain  two  factors  >li(u)  and  a  single  ^st(tt). 
The  other  terms  arc  smaller  by  a  factor  of  e  and  hence  the  Strang  splitting  is  0(ek3) 
accurate.  But  now  suppose  that  on  every  third  step  we  reverse  Si  and  Sa  in  the  Strang 
splitting,  so  that  the  approximate  solution  operator  over  three  timesteps  becomes 

S(U)  m  Sa(k/2)Si{k)Sa(k/2)Si (k/2)S2(k)Sx (k)Sa(k)St (k/2). 

Then  the  0(Jk3)  term  of  the  error  is  simply  twice  the  expression  (2.8)  plus  the  expression 
(2.8)  with  ^i(ti)  and  y!j(u)  interchanged.  The  0(efc3)  terms  then  cancel  leaving  only  the 
0(f3A:3)  terms,  plus  of  course  the  higher  order  terms,  which  are  0(ck*).  Unfortunately  in 
practice  these  higher  order  terms  often  dominate,  especially  when  large  timesteps  k  arc 
used.  Numerical  results  indicate  that  this  modification  has  little  practical  value  except 
when  a  very  fine  mesh  is  used.  This  idea  will  not  be  developed  any  further  here. 
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2.4.  Computing  the  splitting  error  for  quasilinear  systems. 

The  expressions  (2.6)  and  (2.8)  for  the  splitting  errors  look  deceptively  simple. 
Evaluating  them  for  practical  problems  is  actually  quite  a  chore,  mainly  because  of  the 
matrix  derivatives  which  occur.  We  will  now  discuss  the  proper  way  to  evalueate  such 
expressions  and  give  several  examples.  We  are  paticularly  interested  in  the  situation 
where  /l(tt)  =  j4(u)u,. 

We  begin  by  discussing  derivatives  of  matrices.  If  >l(ts)  6  IRrXr  is  a  matrix  valued 
function  of  a  vector  it  €  ntr,  then  its  derivative  A«(tt)  €  IRrXrXr  will  be  a  three-tensor. 
It  is  convenient  to  think  of  this  as  an  array  of  matrices: 


...  f dA  dA  dA] 


(2.12) 


A  tensor  multiplied  by  a  vector  gives  a  matrix.  There  are  two  ways  to  perform  this 
tensor-vector  multiplication  and  it  is  important  to  distinguish  between  them,  since  they 
give  dilTcrcnt  results. 

If  fl  is  the  tensor 

B  =  |fli,fla,...,fl,] 

where  fly  £  IRrXr,  and  if  v  £  lRr,  then  the  first  type  of  multiplication,  denoted  simply 
by  Bv,  is  obtained  by  taking  a  linear  combination  of  the  matrices  fly: 

Bv  =  fli»i  -f  B3V2  +  •  ••  +  Brvr  £  R,Xr 

where  v  =  (i>i, . . . , t»r)T.  The  second  type  of  multiplication  will  be  denoted  by  A(g)v. 
This  product  is  given  by  the  matrix  whose  /th  column  is  the  vector  flyt/: 


Btv\Bav\- ■■  }flrv]  6DtrX< 


fl0V  = 

It  is  easy  to  verify  that  if  w  £  IRr  is  some  other  vector  then 

(B( g)«)w  =  [Bw)v  €  IRr. 


(2.13) 


Roth  forms  of  multiplication  play  a  role  in  diFTerentiation.  Consider  the  vector-valued 
function  f(u)  —  A(u)w,  where  w  is  a  constant  vector  and  u  is  itself  a  function  of  x.  Then 
differentiating  /  with  respect  to  z  gives  the  vector 


s/w = 

=  (A,(u)us)w  6  lRr 


where  /l«( u)  is  the  tensor  (2.12)  and  the  multiplication  is  of  the  first  type.  On  the  other 
hand,  differentiating  with  respect  to  it  gives  the  matrix 


d 


^/(«)  =  /l.(«)®^GlRrXr. 


(2.14) 
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This  can  be  verified  by  directly  computing  the  Jacobian  matrix  corresponding  to 

/(«)  = 

Ei-i  «u(“K 

X;.,  «rj(«K. 

Also  note  that  if  B  is  a  constant  tensor,  then 

a 

—(Butt)  =  i?tt  +  B( g(t». 

Now  let  y?(u)  =  A(u)ux  and  suppose,  for  simplicity,  that  A*(ti)  is  constant,  so  that 
A«v(u)  =  0.  (Otherwise  this  would  be  a  four-tensor.)  We  then  find  that 

^»(u)  =  A„(g)ux  +  Adx 
yivw(u)  =  2A„(jSJ)c?x. 

Example  2.3.  Consider  the  problem  ut  —  ( Af  +  Am(u))ux  with  Af  constant  and  A , 
a  function  of  u  alone.  Take  =  A/u  and  Ai{u)  =  A.(u)ux.  Using  (2.6)  we  can 

compute  the  0(k2)  term  of  the  splitting  error  for  the  first  order  splitting  (1.19): 

%k2(AiuAi{u)  -  Ai«^a(«)) 

=  ^t((iln0ux  +  A,dx)AfUx  —  Afdx(A, t*x)]  ^  jjj 

=  $fca((A.»<g)uxA/iix  +  A,Ajuxx  -  (A/(A„,ttx)ux  +  AfAmuxx)\ 

=  WW  ««A/t»*  —  AfAn  ttx)ux  +  (A,  Af  —  AfA,)uxz\. 

To  obtain  the  last  line  we  have  used  (2.13)  to  rewrite  Aau0uxA/«x  as  A,*A/uxux.  Note 
that  in  order  for  (2.15)  to  be  scro  A/  must  commute  both  with  A,  and  with  A„,. 

As  a  concrete  example,  consider  the  one-dimensional  shallow  water  equations  (1.36) 
with  the  splitting  (1.38).  For  this  system,  the  tensor  A,w  is  given  by 


A,„(u)  = 


1  0 

0  * 

0  1 

i  o 

We  compute  that 
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Similarly,. 

Af(A„ut )t*,  _  |u  ^b  J. 

.  .  .  —  ^o)«*«] 

•  iM/*“  ~  '’'  •*■•  “  2  +  ««..]' 

and  hence 

i**Ota.>li(«)  -  *L*a(«))  =  -  }«.*»  +  )["}]. 

Example  2.4.  The  error  in  the  Strang  splitting  can  be  computed  analogously.  In 
interpreting  the  expression  (2.8)  it  is  important  to  recall  (2.14),  which  indicates  that,  for 
example, 

(A»^a(«))*^i(«)  =  Mi«.<g»fa(u)  +  ^u^8»»Mi(«) 

=  4i»„4i(tt)4a(«)  +  AiuA-2MAi{u). 

Evaluating  (2.8)  for  the  shallow  water  equations  requires  a  tedious  calcuation.  In  view 
of  (1.37),  the  dominant  terms  of  (2.8)  are  the  terms 

$(*u*i(»)Ma(«)  +  K*i»*a(®)Mi(«0  “  Pai.^W)«^iWj- 

This  turns  out  to  be 

~  “*Ul*l  =  0(eVo*5)-  (2.16) 

AH  other  terms  in  the  splitting  error  are  0(c3$}fc3  +  e2$}fc4). 


This  turns  out  to  be 


(2.16) 


2.5.  Efficiency  analysis  for  the  time-split  method  on  hyperbolic  problems. 

The  remainder  of  this  chapter  deals  only  with  hyperbolic  problems  of  the  sort 
described  in  Section  1.4,  although  the  same  type  of  analysis  can  easily  be  applied  to  other 
problems.  For  definiteness  we  also  restrict  our  attention  to  the  Lax-WendrofT  method. 
Other  schemes  can  be  analysed  in  the  same  manner.  In  Section  5.3  a  similar  analysis 
will  be  performed  for  the  Crank- Nicolson  method  on  a  convection-diffusion  problem. 

For  the  constant  coefficient  equation 


ut  =  Avx  —  (A/  +  A,)u, 


(2.17) 


we  wish  to  compare  the  unsplit  method 


LW(A,  It) 


with  the  time-split  method  (1.29)  using 


Q,(k)  =  LW(A.,k). 


{2.18a) 
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For  Q/{kf  2)  we  consider  both 

Qf(k/2)  =  exp( \kA,dK)  (2.18b) 

and 

Qf(k/2)  =  (LW(Af,  k/m)r'*  (2.18c) 

for  some  even  integer  m.  The  split  scheme  defined  by  (2.18a,c)  might  be  used  if  A /  were 
sparse  relative  to  At,  while  (2.18a, b)  would  be  appropriate  for  perturbed  problems  where 
exp(£kA/ds)  is  known  exactly. 

In  each  case  we  assume  that  X  =  k/h  is  fixed  as  Jfc  — ►  0.  In  comparing  the  split 
methods  with  the  unsplit  method,  it  does  not  suffice  to  compare  the  local  truncation 
errors.  For  fixed  k  and  h  the  two  methods  may  take  quite  different  amounts  or  work  to 
implement.  Furthermore  the  optimal  mesh  ratio  may  be  different  for  the  two  schemes. 

Instead  we  compare  the  amount  of  work  required  to  compute  a  solution  with  an 
error  bounded  by  r,  say.  Specifically,  we  consider  the  z-interval  [0, 1]  and  determine  the 
amount  of  work  required  to  compute  solutions  at  time  t  =  1  with  error  no  greater  than 
t.  Strang[49]  takes  an  equivalent  approach  and  comparer  the  accuracy  obtained  with  a 
fixed  amount  of  work.  In  comparing  numerical  results  it  is  convenient  to  take  yet  another 
approach  and  simply  normalize  the  resulting  errors  by  multiplying  by  some  measure  of 
the  work  required  to  obtain  them.  This  will  be  done  in  later  sections. 

For  theoretical  analysis  the  approach  taken  here  seems  to  be  the  most  natural.  It 
determines  the  optimal  mesh  ratio  and  also  provides  (rough)  expressions  for  the  Values 
of  k  and  h  which  must  be  used  to  achieve  a  given  accuracy. 

For  this  analysis  we  will  assume,  as  does  Strang,  that  the  variables  have  been 
normalized  (or  the  norm  appropriately  chosen)  so  that 

p(A)  »  ||A||  =  a 

where  p(A)  is  the  spectral  radius  of  A.  This  means  in  particular  that  ||A3||  a3.  For 

the  splitting  indicated  in  (2.17)  we  suppose  that 

IM/II  » «*  IM.il  «*  (2.19) 

with  the  spectral  radii  again  comparable  to  the  norms  and  f  <  1.  Set  b  =  ca.  Also 
suppose  that  ||w**x||  «  1.  This  is  for  convenience  only,  since  it  removes  one  common 
factor  from  all  of  the  bounds  below. 

Efficiency  of  the  unsplit  method.  We  will  first  analyze  the  unsplit  Lax-Wcndroff 
method  LW(A,k).  Suppose  that  W  is  the  work  required  to  compute  LW[A,k)U J),  at  a 
single  point  xm.  Then  the  work  required  to  advance  the  solution  on  a  unit  z-interval  by 
one  unit  of  time  is  W/kh  —  XW/fc*  if  k  —  \h.  The  truncation  error  for  the  Lax-Wendroff 
method  is  given  by  (1.13), 

ELW(k)u  =  -  $lfc(*2/l3  “  h2A)u„,  +  0(k/i).  (2.20) 

Applying  this  roughly  l/k  times  gets  us  to  time  t  =  1  and 

[LW(A,  k))l'k  =  (exp(fc(A/  +  A.)dx)  +  ELW[k))l'k 

=  exp(A^)  +  ^[ELW{k)  +  0(*'1)]  +  ()(k4). 
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The  error  after  one  unit  of  time  using  the  unsplit  method  is  thus  bounded 

iK(un^  *))'/» -«P(Ad,))«u 

<  j(#‘*IH*ll  +  “*IWI)  +  o(t‘)) 

<  J»V  +  «A*)  +  o(**). 

Since  we  require  an  error  sa  r,  we  set 

i*V+«/x*)  =  T 

giving 


k2  = 


8 r 


o(aa  +  1/X*)* 

Thus  w(t;  X),  the  work  required  to  achieve  a  given  accuracy  r  using  Lax-WcndroiT  with 
mesh  ratio  X,  is  given  by 

.  ..  \W 

«(r;X)  =  -p 

aaW 

=  (Xa  +  l/Xa)— . 

We  have  not  yet  specified  X.  Choosing  X  to  minimise  w(r;  X)  gives  X  =  1/a  and  the 
minimum  work  u/(r)  is 


w(r)  = 


a3W 
3  T 


for  unspiit  Lax-Wendroff. 


(2.21) 


Note  that  the  optimal  mesh  ratio  X  =  1/a  is  also  the  stability  limit  for  this  problem. 
We  can  actually  sec  that  this  is  the  optimal  mesh  ratio  by  looking  only  at  the  error  at 
time  1=1.  Since  this  error  is  bounded  by 


£(*  V  +  h3o)  +  0(k*) 

it  is  clearly  optimal  to  choose  k  and  h  so  that  the  two  terms  k2a3  and  h2a  arc  roughly 
the  same  size  (for  otherwise  wc  could  increase  k  or  h,  and  decrease  the  amount  of  work 
we  do,  without  substantially  increasing  the  error). 

So  far  this  analysis  is  completely  standard  and  our  results  agree  with  those  of 
Slrang[49].  However,  the  same  type  of  analysis,  when  applied  to  time-split  methods 
under  the  assumption  (2.19),  yields  some  illuminating  new  results.  This  will  now  be 
done,  first  for  the  method  (2.18a,b)  and  then  for  (2.18a, c). 

Efficiency  of  the  split  method  (2.1 8a, b)  on  perturbed  problems.  Let  Wt  be 

the  work  required  to  apply  Lax-WendrolT  on  the  slow  scale  and  W™p  the  work  required 
to  compute  exp(M Then  the  work  required  Tor  a  single  step  of  the  time-split 
method  is  WT,M  —  W,  +  2W^xp.  Typically  WT5M  «  W.  The  error  over  one  unit  of 
time  for  the  split  scheme  is  hounded  by 

\\((Q/(k/2)Q.(k)Q/(k/2)fk  -  exp  (Adtj)  «|| 

<  ^ll®»plit(fc)ti  +  lSt(k)u  4-  2/?/(fc/2)u  +  0(fc4)||. 
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For  (2.18b),  E/(k/ 2)  =  0.  The  truncation  error  for  Lax-WcndrofT  on  the  alow  scale  it 
bounded  by 

IUU*H  <  +  M) 

«  1  k*(b*  +  6/X»). 

The  splitting  error  Tor  (2.17)  is  easily  computed  to  be 

£.Piu(*)  =  exp( $*A/d.) exp(kAtd,)  exp(^kAjda)  -  exp(k(Af  +  At)dt) 

=  -  lk*(lA}A.  -  \AfA,Af  +  {A.A)  (2.22) 

-  \A\Aj  +  A.A,A.  -  \A,A\)d\  +  0{k*) 

so  that 

l|£.piit(*H|  <  i k'ia'b  +  ab*)  « 

although  it  may  be  much  smaller  for  some  problems.  Since  our  results  depend  very  much 
on  the  site  of  this  error,  we  will  suppose  for  now  that 

P'.Piu(*HI  < 

for  some  a,  so  that 


i||^.piu(Jk)«  +  E.(k)«\ |  <  $lfc2(o  +  6*  +  fr/X2). 
In  order  to  obtain  accuracy  r  we  must  take 

.2 _  6t 

o  +  b3  +  6/X2 


w(t;X)  =.XVFTaM/fc2 

=  X(«t  +  63  +  6/X2)^~— — .  (2.23) 

OT 

The  optimal  stepsise  ratio  X  now  depends  on  the  site  of  the  splitting  error  and  is  given 

*>y 

x  -  ‘22,> 

so  that 


w(t)  =  \/6(<7  +  &3)  — ^ —  for  the  time  split  method  (2.18a, b). 

•jT 

If  a  <  6s  (e.g.,  when  Af  and  A,  commute),  then  (2.24)  gives  X  1/6  and 


w(t) 


63H'T5M 

3t 


(2.25) 
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Table  2.1 


Reduction  in  work  over  (2.21 )  obtained  by  using  the  time-split  method 
(2.18a, b)  on  (2.17).  The  results  depend  on  the  site  of  the  splitting  error. 


case 

£*piu(fc) 

optimal  X 

reduction  in  work 

general 

o/t* 

/  " 

V  0  -f  e*a* 

>/«(*  +  c3*s) 

best 

0 

1 

ea 

ea 

typical 

ea3*3 

1 

e 

a 

When  W™*  *=s  W  this  is  better  than  (2.21)  by  a  factor  of  ea,  meaning  greatly  improved 
efficiency.  Note  that  when  o  —  0  the  only  error  incurred  is  the  error  in  using  Lax- 
WcndrofT  on  the  slow  scale.  From  our  previous  discussion  of  Lax-Wcndroff  it  is  clear  why 
X  =  l/b  is  optimal  in  this  case. 

On  the  other  hand,  if  the  splitting  error  i?  as  bad  as  (2.22)  indicates,  then  o  —  o26 
and  X  1/a  in  (2.24),  giving 


w(t)  = 


obWTgM 

3t 


This  is  still  an  improvement  over  (2.21),  although  now  by  only  a  factor  of  c.  Note  that 
now  X  is  chosen  appropriate  to  the  fast  scale,  even  though  the  fast  part  of  the  problem 
is  solved  exactly.  This  is  necessary  because  of  the  splitting  error.  Indeed,  if  we  try  to 
use  X  =  l/b  when  o  =  aah,  we  obtain  no  improvement  over  (2.21).  For  this  reason  it 
is  advisable  to  always  use  small  timesteps  with  the  time-split  method  (2-J8a,b)  unless 
E,p[\i(k)  is  known  to  be  very  small,  in  which  case  even  greater  efficiency  is  achieved  by 
using  larger  timesteps. 

These  results  arc  summarized  in  Table  2.1. 


Efficiency  of  the  split  method  (2.18a,c)  with  sparse  Aj.  When  Lax-WendrofT 
is  used  Tor  both  operators,  the  work  for  a  single  step  of  the  lime-split  method  is  given 
by  VFT*M  =  W,  +  iriWf,  where  Wj  is  the  work  required  to  apply  Lax-WcndrolT  on  the 
fast  scale.  We  arc  assuming  that  Wf  <  Wm  s=3  W.  Suppose  that  Wf  =  7 IF  for  some 
7  <£  1.  In  this  case,  the  best  we  can  hope  for  is  to  decrease  the  required  work  by  a  factor 
of  7.  Wc  will  see  that  in  general  we  can  reduce  the  work  by  a  factor  of  roughly  7  -f  \/e 
by  choosing  the  mesh  ratio  appropriately.  When  the  splitting  error  is  negligible,  wc  can 
improve  this  to  7  +  f. 

Ikcause  wc  arc  still  free  to  choose  m  in  (2.18c),  the  mesh  ratios  we  use  on  the  fast 
and  slow  scales  arc  essentially  independent  for  this  problem.  The  local  truncation  error 
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for  (2.18e)  U 

Ef(k/2)  =  {IAY(Af,k/m)r't -exp^kAfd,) 

=  («P<;M/*«)  -  Ht&Af  -  mh2Af)dl)m/t  -  np(^kAfda) 

=  -i{£'A,/-ih*A/)dl+o(k*). 

The  optimal  value  of  m  is  that  which  makes  I t*a*/ma  s=s  h2,  or  m«  Xa  where  X  —  kfh 
is  the  mesh  ratio  for  the  slow  scale.  The  optimal  mesh  ratio  on  the  fast  scales  is  thus 
k/mh  =  1/a  regardless  of  X. 

Using  this  value  of  m,  we  compute  the  following  bound  for  the  error  at  time  t  =  1, 
using  (2.1), 


i||/?.p»t(*)«  +  E.(k)u  +  2Ef(k/2)u\\  <  J**(o  +  6s  +  b/\ 2  +  a2/m2 

fc*(cr  +  6s  +  2a/Xa). 

We  then  obtain 


+  o/Xa)  +  0(ks) 
(2.26) 


w(t;  X)  X(o  +  h3  +  2 a/Xa)^*  ^  for  (2.4a, c). 

or 


(2.27) 


The  optimal  X  is  most  easily  determined  by  requiring  that  the  terms  in  the  error  (2.26) 
balance.  This  gives 


(2.28) 


Again  we  will  consider  the  best  and  worst  cases,  a  —  0  and  a  =  a26.  When  the  splitting 
error  is  negligible,  (2.24)  gives 


(2.29) 


In  this  case  the  optimal  mesh  ratio  appears  to  be  larger  than  the  optimal  mesh  ratio 
for  the  slow  problem  alone  (which  would  be  1/ca).  This  counterintuitive  result  is  due  to 
the  fact  that  otherwise  the  error  on  the  fast  scale  dominates  the  error  on  the  slow  scale, 
lly  taking  larger  timesteps  on  the  slow  scale  we  decrease  the  work  without  increasing 
the  error,  or  so  the  efficiency  analysis  tells  us.  Unfortunately,  the  mesh  ratio  (2.29)  is 
larger  than  the  stability  bound  for  LW(A,,k),  which  is  1/ca,  and  so  this  cannot  be  used 
in  practice.  The  best  wc  can  do  is  to  take 

X«s  ~ 
eo 


with  corresponding  work 


w(r)  = 


2 e2a3  Wt  +  c-lW/ 
co  6r 


=  2a3j»L-t*> 

8  r 


(<  +  7) 


a2W 

3t 
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case 

optimal  X 

reduction  in  work 

general 

ah3 

minf 1  J  2a  ] 

_ ( _  .  Iff  +  e3a3\  , 

yea'i  e  +  t3a3 J 

max  1  t,  y  1  +  7 

best 

0 

1 

£  +  7 

to 

typical 

to3*3 

1 

+7 

y/ea 

which  is  better  than  (2.21)  by  a  Factor  of  £(t  -f  7). 

In  the  more  typical  situation,  when  a  =  a2b,  (2.28)  becomes 


X  « 


with  corresponding  work 


w(t)  = 


3eo3 

y/la 


W,  + 

6  T 


_a,y/tW'  +  Wf 
2  T 

«(>A  +  7)-^r. 

Wc  thus  see  that  if  yfl  <  7,  an  increase  in  efficiency  by  the  best  possible  factor  of 
roughly  7  is  always  possible.  These  results  arc  summarized  in  Table  2.2. 


2.6.  Phase  errors. 

When  solving  differential  equations  with  wave-like  solutions,  it  is  frequently  desirable 
to  compute  the  phase  errors  of  the  finite  difference  scheme  employed.  Comparing  the 
phase  errors  for  the  lime-split  method  with  those  for  the  unsplit  method  provides  some 
further  insight  into  the  results  of  Section  2.5. 

Consider  again  the  constant  coefficient  problem  ut  =  Aux  and  denote  the  eigenvalues 
and  eigenvectors  of  A  by  py  and  «y,  respectively, 

AUj  =  /ty«y,  j  =  1,2,  ...,r 
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with  the  py  ordered  as  in  (1.26).  As  usual  we  suppose  that  ||A||  r=s  p(A)  =  pr.  If  we  take 
as  initial  conditions  a  single  mode 


tt(z,  0)  =  e'^tty 


(2.30) 


for  some  j  and  some  wavenumber  then  the  true  solution  at  time  t  is  simply 


The  wave  thus  propagates  with  a  phase  speed  py. 

Now  suppose  we  apply  a  single  step  of  unsplit  Lax-Wcndroff  to  u(xt0).  By  (1.13)  we 
obtain 


LW(A,  k)u(x,0)  =  u(x,  k)  -  £jfc(fc2A3  -  h2A) ttx«(z,0)  +  0{k*) 

=  -  |4(*X  -  fcVy)(*03K  +  0(k *) 

=  exp{i£(z  +  k(pj  +  £fc2(p3  -  MyA2)£2))}fiy  +  0(*4). 

The  phase  speed  of  the  numerical  wave  is 

My  +  +  0(fc3). 

The  optimal  mesh  ratio  for  Lax-Wendroff  is  X  »  l/)pr|.  In  practice,  of  course,  one  never 
has  exactly  the  optimal  mesh  ratio,  so  we  suppose  only  that  X  =  1/p  with  p  >  |pr|.  We 
then  find  that  the  error  in  the  phase  speed  for  the  j th  eigenvector  with  wavenumber  (  is 

phase  speed  error  =  g £2(p3  —  PyP2)f2  +  0(£3). 

For  comparison  purposes  we  again  wish  to  normalize  by  some  measure  of  the  work 
required  to  compute  the  solution.  We  define  the  normalized  phase  speed  error  ^y(£)  as 

=  (phase  speed  error )/kh. 

For  the  unsplit  method  we  have 


*>(*)  =&My-MyM2)*2 +  <?(*)■ 

If  py  =  p  exactly  then  there  is  no  error  in  this  mode  of  the  computed  solution.  In 
general,  however,  the  error  is  roughly 

MO  «  -  b^r?.  (2.31) 

Now  consider  the  split  method  (2.18a, b)  where  exp '  \kAjdT)  is  known  exactly  and 
suppose  to  begin  with  that  there  is  no  splitting  error  for  the  splitting  A  —  Af  +  Aa. 
Then  the  matrices  are  simultaneously  diagonalizablc  and  so  the  tty  arc  also  eigenvectors 
of  Af  and  A,.  We  then  have 

A,«y  —  Pay  «y 
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with  <  (\pT\-  The  optimal  mesh  ratio  as  found  in  Section  2.5  is  then  X  «  1/e/t. 
When  applying  (2.18a, b),  the  only  error  is  the  Lax-Wendroff  error  on  the  slow  scale,  so 
after  applying  one  step  of  the  split  operator  Q(k)  we  obtain 

Q(*)tt(z,0)  =  k)  -  bk(kaA]  -  h2As)usxx[x,Q)  +  0(k*). 

Proceeding  exactly  as  before  we  find  that  the  normalized  error  is 

MO  =  -  fV«iM2)£2  +  0{k) 

»  -  i<P.j/*r{2 

This  is  always  better  than  (2.31).  Just  how  much  better  it  is  will  depend  on  the  velocity 
of  the  mode  (2.30).  For  slow  waves,  those  for  which  |/i,|  <  c|pr|,  say,  we  have  \paA  pa 
l/Xjj  and  so  (2.32)  is  better  than  (2.31)  by  roughly  a  factor  of  t.  For  fast  waves,  on  the 
other  hand,  for  which  |/iy|  ps  |/tr|,  (2.32)  is  better  than  (2.31)  by  a  factor  of  t2.  The 
improvement  in  phase  errors  is  thus  more  dramatic  for  fast  waves  than  for  slow  waves. 
This  is  to  be  expected  since  it  is  the  fast  subproblem  which  is  being  solved  exactly. 

How  do  these  results  fit  in  with  the  results  of  Section  2.5?  There  we  saw  that  for  the 
method  (2.18a,b)  with  no  splitting  error,  the  work  required  to  obtain  a  given  accuracy 
should  be  reduced  by  e2,  or,  equivalently,  that  the  normalized  error  should  be  reduced 
by  c2.  Yet  here  it  seems  that  the  error  in  slow  waves  is  reduced  only  by  t.  This  apparent 
contradiction  is  resolved  by  reexamining  (2.31).  This  shows  that  for  the  unsplit  method 
phase  errors  in  slow  waves  arc  already  smaller  by  a  factor  of  e  than  those  in  fast  waves. 
Hence  with  the  unsplit  method  errors  in  the  fast  waves  dominate,  and  reducing  those 
errors  by  c2  (and  errors  in  the  slow  waves  by  c)  causes  the  overall  global  error  to  decrease 
by  t2. 

This  has  an  important  consequence  which  was  not  directly  apparent  from  the  analysis 
of  Section  2.5.  For  problems  in  which  fast  waves  arc  absent  from  the  solutions  of  interest, 
and  only  slow  waves  arc  present,  the  use  of  the  time-split  method  can  be  expected  to 
decrease  the  normalized  errors,  and  hence  improve  the  efficiency,  by  at  most  a  Tactor  of 
t,  even  in  the  absence  of  splitting' errors. 

Now  suppose  that  the  splitting  error  is  nonzero.  For  tire  constant  coefficient  system 
this  means  that  A /  and  A,  do  not  commute  and  the  eigenvectors  tty  of  A  are  no  longer 
eigenvectors  of  Af  and  As.  Because  of  this  initial  conditions  consisting  of  a  single  mode 
(2.30)  no  longer  lead  to  a  single-mode  solution  and  we  are  not  able  to  consider  each  mode 
separately. 

Instead  we  take  more  general  initial  conditions 

u(x,0)  -  c'*xu 

where 

r 

U  =  ^  (2.33) 

m==  1 

and  look  at  phase  errors  in  the  j th  mode.  We  assume  that  the  ttm  arc  order  unity  and 
for  convenience  take  otj  =  l. 


The  truncation  error  for  the  split  method  is  now  the  sum  of  the  truncation  error  for 
Lax-Wendroff  on  A,  and  the  splitting  error  (2.22).  So 

Q{k)u(x,  0)  =  u(z,  *)  -  $  kBuxxx  +  •  •  •  (2.34) 

whcFB 

B  =  [k*A]  -  h*A.)  +  k2[\A2fA,  -  %AfA,Af  +  \A,A) 

-  \A\Af  +AaAfA.  -  bAfAl). 

Assuming  as  usual  that  ||A,||  <  c||A/||  «  tfiT  and  using  the  optimal  mesh  ratio  k/h 
l/|/ir|  gives  a  rough  bound  on  B : 

||fi||  <  2 k2e£  +  0[k2 e2tf).  (2-35) 

Now  we  must  make  an  additional  assumption  on  the  matrix  A,  namely  that  the 
eigenvectors  of  A  are  well-conditioned.  If  X  is  the  matrix  of  eigenvectors  um  then  we 
assu  me 

||X||  j|X-'||  =  0(l). 

This  means  that  we  can  expand  Bit  as 

r 

B»=  Y  /3mum 

m=l 

with  \/3m\  of  the  same  order  as  ||£?||.  This  is  because  ft  =  X~lBXa  and  so  ||/?||  < 
||/?||  ||X||  ||X-1||  ||q||  «  ||/j||.  Using  (2.33)  in  (2.34)  then  yields 

Q(fc)u(z,0)  =  Y  ame‘C(*+Mmfc)t «m~i*(*€)3  Y  Pmum  +  0{k4) 

m=sl  m=l 

r 

=  Y  a’n  +  6  ^2^m/«m)]}«m  +  0{k*). 

m  =  l 

Using  (2.35)  we  can  compute  the  normalized  phase  speed  error  in  the  jth  mode: 

*»(€)  «  Mt*- 

Note  that  unlike  the  previous  cases  the  phase  error  here  is  the  same  for  fast  waves  and 
slow  waves.  Comparing  this  with  (2,31)  shows  that  for  fast  waves  [ft}-  s=»  fir)  the  error  is 
reduced  by  e  while  for  slow  waves  (fij  <  cfiT)  the  error  is  not  reduced  at  all.  This  indicates 
that  when  computing  a  solution  containing  only  slow  waves,  the  time-split  method  with 
splitting  errors  may  be  no  more  efficient  than  the  unsplit  method. 


2.7.  Block  triangular  systems. 


Since  the  efficiency  or  the  split  scheme  is  limited  primarily  by  the  splitting  error, 
it  is  interesting  to  investigate  how  this  error  depends  on  the  coupling  between  fast  and 
slow  scales  in  a  simple  model  system.  Consider  the  block  triangular  system  with 


An 

An 


and  the  splitting 

0  An 
0  Aa  a 

and  suppose  that  ||/lyy||  1  and  that  ||j4i2||  =  a  <  1.  For  variety  we  have  chosen  a 

problem  in  which  Aj  — ►  oo  as  £  — ►  0  rather  than  A,  —*  0.  The  theory  developed  in  the 
previous  section  applies  equally  well  in  this  situation. 

Here  A 12  is  the  coupling  between  fast  and  slow  scales.  If  A\a  =  0,  the  problem  is 
uncoupled  and  Eap i;t(fc)  =  0.  In  general,  from  (2.22), 


Af  = 


Mu  0 

0  0 


A,  = 


7£*plit(fc)  —  g 


^Mn(Mu^ia  ~  ZAiaAn)  g3  +  o(Jfc4). 


Thus  ||^,piit(l:)u||  «=s  «fc3/24c2.  The  efficiency  of  the  splitting  depends  on  the  size  of  a. 
In  the  notation  used  above,  we  have 

a=l/c,  6=1,  a  —  \aa3b. 

For  unsplit  Lax-Wendroff,  (2.21)  gives 


wiT)  =  (2-36) 

The  time-split  method  (2.l8a,b)  is  always  more  efficient  if  we  choose 

X  ss  (1  +  %aa2b)~l/2. 

For  example,  if  a  1  we  should  use  X  «  2/a  =  2c  in  order  to  reduce  (2.3ft)  by  a  factor 
of  c.  The  maximum  efficiency  indicated  in  (2.25)  is  achievable  only  if  a  <  ca,  in  which 
case  taking  X  =  1  reduces  (2.36)  by  a  factor  or  cs. 

2.8.  Reducing  the  splitting  error. 

For  block  triangular  systems  in  which  /1|2  is  not  small,  it  is  possible  to  reduce  the 
coupling  through  a  change  of  variables  so  that  the  optimal  efficiency  can  be  achieved.  A 
change  or  variables  amounts  to  replacing  u  by  ti  =  13u  for  some  nonsingular  matrix  B. 
The  system  ut  =  Aux  then  becomes  ut  —  BAU~xux.  Clearly,  if  71  is  chosen  to  be  the 
eigenvector  matrix  of  A  then  the  problem  completely  decouples  into  independent  scalar 
equations.  We  arc  seeking  something  less  expensive  which  only  decouples  the  fast  and 
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slow  scales.  Thus  we  want  a  (well-conditioned)  matrix  B  such  that 


BAB-1  = 


(2.37) 


with  ||Cn||  *=s  ||Cja||  »  1*  In  the  block  triangular  case,  it  suffices  to  consider  B  of  the 
form 

T  f  r%  "1 

B-1  = 


Then 


BAB 


>  _  \l  Bial  n-i  I  -Bial 

lo  /  ]’  [o  /  J- 

"=[!o" 


—  jAuBia  +  An  +  B12A32 

Am 


and  so  Bn  should  be  chosen  to  solve 


-AnBn  —  B12A22  =  An 


(2.38) 


in  order  to  completely  decouple  the  fast  and  slow  scales. 

In  the  present  context  solving  for  Bn  from  (2.38)  is  not  worthwhile.  In  order  to 
achieve  optimal  efficiency  we  need  only  reduce  the  coupling  by  one  or  two  factors  of  c. 
Further  reductions  do  not  gain  anything  once  the  Lax-Wendroff  errors  dominate.  This 
suggests  taking 

Bn  ~  {Aii  An  (2.39) 


so  that 


where 


BAB'1  =  \*An  Ara 
i  0  /I22 . 

a\%  ~  tAi1  AnA22- 


We  now  have  U-dj^U  ca  provided  ||/17|,||  **  1.  The  coupling  is  thus  reduced  by  a 
factor  of  e  through  the  use  of  a  very  simple  change  of  variables.  This  process  can  be 
repeated  to  obtain  additional  factors  or  c.  This  change  of  variables  has  been  suggested 
by  Kreiss[32]  in  a  similar  context. 

For  full  systems  of  the  form 


An 

[  .A21  An 

we  can  obtain  a  similar  reduction  in  the  sir.e  of  both  off-diagonal  blocks  and  again  reduce 
the  splitting  error  by  several  orders  of  magnitude.  In  this  case  we  consider  li  of  the  form 


H=J'  K]\l  °1  _  [/  +  KIj  K 

LO  /JU  /]  i  I 


It  is  easy  to  verify  that  the  lower  corner  of  A  is  annihilated  by  taking  L  to  satisfy 

—IjA  1 1  —  A22L  -  LAnL  +  A21  —  0. 
t 
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The  matrix  K  can  then  be  chosen  as  before  to  remove  the  remaining  upper  corner.  This 
results  in  a  system  of  the  form  (2.37).  This  particular  transformation  is  discussed  more 
completely  by  O’Malley  and  Anderson[44j.  Again,  however,  we  are  not  interested  here  in 
completely  annihilating  the  corners,  but  rather  in  reducing  them  by  a  factor  of  c.  This 
is  easily  accomplished  by  taking 

K  =  eATtlAi2 
L  —  -  (A2tAif . 

Example  2.4.  This  problem  is  designed  to  illustrate  the  effects  of  the  splitting  error 
and  the  use  of  the  change  of  variables  (2.30).  Consider 

[“  ^  =  l0°  j  [“  for  0  <  x  <  1,  t  >  0  (2.40) 

with  initial  conditions 

t»(x,0)  =  v{x,0)  =  c-,0<,(*-»/2)a 
and  periodic  boundary  conditions 

u(0,t)  =  u(t,t),  t  >  0,  /  =  1,2, 

v(O,0  =  v(i,<),  <>0,  y  =  l,2. 

Figure  2.1a  shows  the  results  after  236  time  steps  using  Lax-Wendroff  with  h  =  1/50  and 
k  =  A/10  on  the  unsplit  problem.  Figure  2.1b  shows  the  results  based  on  the  splitting 


We  used  k  =  h  —  1/50  with 

Qa(k)  =  LW(A„  k),  Qf(k/ 2)  =  (LW(Af,  A/10))5. 

In  this  case  Ea(k)  =  Ef(k/ 2)  =  0  by  a  judicious  choice  of  k/h  and  m.  The  second 
component  v  is  computed  exactly  and  the  errors  in  u  are  due  entirely  to  the  splitting 
error. 

If  the  change  of  variables  suggested  in  (2.39)  is  applied  twice  to  (2.40)  with  e  s  0.1, 
we  obtain  the  new  variable 

u  =  u  —  (e  +  ca)v  =  u  —  O.llv  (2.41) 

and  (2.40)  becomes 


If  we  solve  this  system  with  the  same  split  scheme  as  before  and  then  transform  back  to 
the  original  variables  by  «  =  u  +  O.llv,  the  errors  in  u  are  reduced  to  O(10~3)  as  seen 
in  Figure  2.1c. 
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Fig.  2.1.  True  ( dashed  Une)  and  computed  solutions  at  t  =  4.72  for  Example 
2.1.  The  first  component,  u,  is  on  the  left  and  the  second  component,  v, 
is  on  the  right.  The  schemes  used  are:  (a)  unsplit  Lax-Wendroff,  (b)  the 
time-split  method  (2.18a, b),  and  (c)  the  time-split  method  with  the  change 
of  variables  (2.41). 


2.9.  The  shallow  water  equations. 

In  this  section  the  efficiency  analysis  of  Section  2.5  is  applied  to  the  one-dimensional 
shallow  water  equations  (1.36).  We  will  use  the  splitting  (1.38)  and  assume  that  the 
condition  (1.37)  holds.  In  Section  2.4  we  computed  the  splitting  error  for  this  system. 
In  general  this  error  is  nonncglible.  Since  all  of  the  waves  in  the  solution  to  the  original 
problem  are  fast  waves,  the  analysis  of  Sections  2.5  and  2.6  leads  us  to  expect  the  time- 
split  method  to  be  more  efficient  than  the  unsplit  method  by  a  factor  of  t. 

Taking  the  mesh  ratio  as  in  (1.39),  the  time-split  method  (2.18a,b)  becomes 


U 


im 


tC+1 

*n+t 


=  MtC-r  +  Unm+p  +  +p] 

=  ^-r-^+p  +  ^-P  +  ^+rl 


=  LW(A.,k) 


=  flic,  +  tC+p  +  C-P  -  <+p) 
=  r;*-p  -  <c+r + + <+P\- 


(2.42) 


Since  ut  —  A„ux  is  a  quasilincar  problem,  an  appropriate  generalization  of  the  Lax- 
WcndrolT  operator  must  be  used  for  LW(Aa,k).  We  have  used  MacCormack’s  method 
(sec  [41]). 

We  wish  to  compare  the  efficiency  of  the  split  method  with  that  or  the  unsplit 
method.  For  convenience  in  checking  our  predictions  against  experimental  results,  we 
choose  to  compare  the  error  at  a  fixed  time  normalized  by  the  amount  of  work  required 
to  compute  it  (rather  than  the  amount  of  work  required  to  compute  a  solution  with  a 
given  error).  Since  the  split  and  unsplit  methods  take  roughly  the  same  amount  of  work 
per  grid  point  per  timestep,  it  suffices  to  normalize  the  errors  at  a  fixed  time  by  dividing 
by  kh,  as  we  did  to  normalize  the  phase  errors  in  Section  2.6. 

We  first  consider  the  unsplit  MacCormack’s  method  applied  to  (1.36).  Since  A  «  Aj 
with  small,  slowly  varying  pclurbations,  the  errors  in  applying  MacCormack’s  method 
on  A  arc  roughly  the  same  as  those  in  applying  Lax-Wendroff  on  the  constant  coefficient 
matrix  A/.  We  can  thus  use  the  results  of  Section  2.5  directly  to  analyze  the  efficiency 
of  the  unsplit  method. 

Since  p[A)  p(Af )  =  $0/2,  the  optimal  mesh  ratio  is  X  «  2/</>0.  The  error  at  time 
t  =  1  is  bounded  using  the  truncation  error  (2.20)  by 

i||JS"*r(*)fi||  «  l(k2||yts||  +  A2||4||)  ||3i«||.  (2.43) 


For  smooth  solutions  we  can  assume  that  j|uxxx||  =  0((<j> o).  Then  taking  X  =  O(l/^o), 
we  (ind  the  normalized  error  by  dividing  (2.43)  by  kh: 

normalized  error  =  O(c^jj)  for  the  unsplit  method.  (2-44) 


Now  to  analyze  the  split  method.  The  splitting  error  (2.22)  is  in  general 
for  this  problem.  The  results  of  Section  2.5  indicate  that  for  the  time-split  method  with 
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Ej  =  0  and  Etp\n(k)  non  negligible,  we  should  again  take  X  =  O(l/0o)  and  hence  the 
optimal  p  in  (1.30)  should  be  some  small  integer,  independent  of  both  e  and  0o.  Numerical 
experiments  confirm  this  prediction  (sec  Example  2.5  below)  and  in  fact  p  =  3  or  4  seems 
to  be  optimal  over  a  wide  range  of  values  of  t  and  ho¬ 
using  this  optimal  value  of  X,  the  normalized  error  should,  in  theory,  be  reduced  by 
a  factor  of  e  over  (2.44),  i.e., 

normalised  error  =  O(£30o)  for  the  split  method.  (2.45) 

This  is  also  confirmed  in  the  following  example. 

Example  2.5.  Consider  the  shallow  water  equations  (1.36)  on  0  <  x  <  1  with  initial 
conditions 

tt(z,0)  =  e<f>  qcos(2xx) 

0(x,  0)  =  0o(l  +  e  sin(2irz)) 
and  periodic  boundary  conditions 


u(0  ,t)  =  u(l,t) 

0(0,1)  =  0(1,<). 

We  first  compare  the  error  obtained  at  a  fixed  time  using  various  values  of  p  in  the 
time-split  method  (2.42).  Figure  2.2  shows  the  normalized  errors  as  a  function  cf  p  for 
0o  =  1  and  c  =  10~3,  10~3, 10-4  with  h  —  1/50.  Other  values  of  0O,  £,  and  h  have  also 
been  tested  and  lead  to  graphs  which  are  qualitatively  very  similar  to  Figure  2.2.  In  all 
cases  p  =  3  or  4  is  optimal. 

We  can  also  compare  the  error  in  the  split  method  with  that  of  the  unsplit  method 
using  the  optimal  values  of  X  for  each.  For  the  split  method  we  take  p  —  3  (corresponding 
to  X  =  12/0o )  and  for  the  unsplit  method  we  use  X  =  l/0o-  Figure  2.3  shows  the  results 
for  0o  =  1.  We  sec  the  normalized  error'  plotted  as  a  function  of  e.  This  confirms 
the  prediction  that  using  the  split  method  reduces  the  normalized  error  by  a  factor  of 
c.  More  significantly,  it  show  that  even  for  i.  irly  large  (i.e.  realistic)  values  of  c  the 
time-split  method  is  superior.  For  example,  at  c  =  0.1  the  errors  are  reduced  by  a  factor 
of  roughly  100. 

Simple  waves.  The  splitting  error  for  the  quasilinear  problem  (1.36)  with  the 
splitting  (1.38)  depends  on  u  and  0  and  the  relation  between  them.  In  general  it  io 
nonnegligiblc  but  for  certain  special  solutions,  namely  simple  waves,  the  splitting  error 
is  identically  zero. 

The  equations  (1.36)  can  be  written  in  characteristic  form  as 

(u  +  0),  =  -($0  +  u)(u  +  0), 

(w  -  0)t  =  ( £0  -  «)(«  -  0)*- 

The  Riemann  invariants  u  +  0  and  u  —  0  arc  each  constant  along  characteristic  curves 
in  x-t  space  defined  by  the  ordinary  differential  equations 

^  =  l<f>(x,  t)  J,  u{x,t) 
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FIG.  2.2.  Normalized  errors  in  the  shallow  water  equations  of  Example  2.5 
as  a  function  of  the  parameter  p  occurring  in  the  mesh  ratio  (1.39).  In  all 
cases  t  =  0.98,  <t>o  =  1  and  h  —  1/50. 


Fig.  2.3.  Normalized  errors  in  the  shallow  water  equations  of  Example  2.5  as 
a  function  of  e  for  the  unsplit  method  with  X  =  \/4>o  and  the  split  method 
with  X  =  12/^o-  In  the  computations  shown  here  t  =  0.96,  0O  =  l  and 


and 


jjjj  =  -(s 0) 

respectively. 

A  solution  for  which  one  ofjthe  invariants  is  in  fact  constant  for  all  z  and  t  is  called 
a  simple  wave.  For  simple  waves  the  splitting  error  is  identically  zero.  This  is  most  easily 
seen  by  changing  variables.  Set 


p(x,t)  =  u{x,t)  +  <l>(x,t), 
o{x,t)  =  u {x,t)-<f>{x,t). 


(2.47) 


The  equation  (1.36)  becomes 


_  _  i 
—  ? 


3  p  +  j  0 

0  p  +  3a 


(2.48) 


The  matrix  occurring  here  is  SAS  1  where 


S  = 


1  1 
1  -1 


Applying  the  same  similarity  transformation  to  A/  and  A,  leads  to  the  following  splitting 
of  (2.48): 


SAfS =  $ 


®  1,  SA.S-1  =-\\*p  +  24>0  .  ,°  .  , 

0  ^oj  *  0  P  +  3<r  +  - 


(2.49) 


Since  we  have  applied  a  constant  similarity  transformation,  it  is  easily  verified  that 
the  splitting  errors  corresponding  to  the  splitting  (1.38)  and  (2.49)  are  also  related  by 
the  same  similarity  transformation.  Thus  it  suffices  to  show  that  for  simple  waves  the 
splitting  error  in  (2.49)  is  zero.  This  is  easy  to  do,  as  we  will  sec  momentarily. 

Wc  note  in  passing  that  solutions  to  the  shallow  water  equations  can  be  computed 
directly  in  terms  of  p  and  a  using  the  splitting  (2.19).  With  ft  and  S  denoting  approxima¬ 
tions  to  p  and  a,  the  time-split  method  (2.42)  then  becomes 


=  /C-. 

5^  =  Sm+l 

n 


=  W(A.,k) 


/p"+ 1  _  if** 

nm  —  'cm-l 
t»n+l 


ft 

is. 


sz+l  =  s 


m+1  * 


(2.50) 


This  form  will  prove  particularly  convenient  when  specifying  boundary  conditions  for  the 
intermediate  solutions,  as  we  will  sec  in  Section  4.5. 

Suppose  now  that  we  are  computing  simple  waves  and  that  one  of  the  invariants  p 
or  a  is  constant,  say  a  —  —fa.  Then  clearly  S JJ,  =  —<t>o  in  (2.50)  and  so  the  second 
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component  of  the  splitting  error  is  zero.  The  equation  for  p  is  pt  =  —  }$p+o)px  =  A(p) 
and  it  remains  only  to  show  that  there  is  no  error  in  using  the  splitting  A\(p)  =  —  \<i>opx, 
Ai(p)  =  —  ^3p  +  c  —  2<f>o)px.  Since  a  and  are  constant,  this  is  essentially  the  problem 
of  Example  2.1(c)  and  so  the  splitting  error  is  zero.  Note  that  the  expression  (2.16)  is 
consistent  with  this,  since  for  simple  waves  uz  =  <f>x  and  uxx  =  tj>xx. 

It  follows  that  the  optimal  mesh  ratio  for  computing  simple  waves  is  0(1  /e^o)  leading 
to  normalized  errors  which  are  reduced  by  a  factor  of  r3  over  (2.44): 

normalized  error  =  0(e3^)  for  the  split  method  on  simple  waves. 

These  predictions  are  also  confirmed  by  numerical  experiments,  as  the  following  example 
shows. 

Example  2.6.  Consider  the  shallow  water  equations  (1.36)  on  0  <  x  <  1  with  initial 
conditions 

u(x,0)  =  t^osin(2xx) 

^(x,0)  =  <£o(l  +  csin(2xx)) 

and  periodic  boundary  conditions.  Since  u  —  ^  is  constant,  the  solution  is  a  simple  wave. 

We  again  compare  the  normalized  errors  at  a  fixed  time  using  various  values  of  p  in 
the  time-split  method  (2.42).  We  expect  p  =  0(1/ e)  to  be  optimal.  In  order  to  test  this 
theory  when  i  is  small  we  must  run  the  computations  out  to  large  times,  t  =  0(l/c).  For 
each  value  of  c  we  will  compare  the  normalized  error  at  t  =  0.96/(100c),  using  values  of 
p  <  12/(100f).  This  is  roughly  the  stability  limit  of  the  method.  (In  Section  3.5  it  will 
be  shown  that  the  stability  limit  is  k/h  <  l/(2<^o)  which  corresponds  to  p  <  l/(8e).) 
Since  the  stability  limit  is  smaller  than  the  optimal  p  predicted  by  the  theory,  we  expect 
the  normalized  errors  to  be  monitonically  decreasing  up  to  the  stability  limit.  This  is 
confirmed  in  Figure  2.4. 

The  theory  also  predicts  that  the  residing  normalized  errors  at  a  fixed  time  should 
be  0(c3<f>l)  at  the  optimal  p,  and  hence  that  the  errors  at  time  0(l/e)  should  be  0(e3^Q). 
This  is  also  confirmed  by  Figure  2.4. 
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Fig.  2.4.  Normalised  errors  in  a  simple-wave  solution  to  the  shallow  water 
equations  of  Example  2.6  as  a  Function  of  the  parameter  p  occuring  in  the 
mesh  ratio  (1.39).  In  all  cases  4>o  —  1  and  h  =  1/50  while  t  —  0.06/(l00c). 


$.  Cauchy  stability 


3.1.  Introduction  to  stability  theory. 

In  this  chapter  we  investigate  the  stability  of  the  one-dimensional  time-split  method 
when  applied  to  a  constant  coefficient  hyperbolic  problem  on  the  entire  real  line,  —  oo  < 
i  <  oo,  or  on  a  finite  interval  with  periodic  boundary  conditions. 

We  will  first  engage  briefly  in  a  general  discussion  of  Cauchy  stability  for  a  marching 
scheme  of  the  form 

t/n+l  =  Q[k)Un  (3.1) 

appled  to  a  constant  coefficient  problem.  More  details  can  be  found  in  Richtmycr  & 
Morton(46]  or  Thom£c(51|.  We  use  a  standard  definition  or  stability,  which  can  be 
written  in  several  equivalent  forms.  We  begin  with  the  most  natural  of  these. 

STABILITY  Definition  3.1.  The  operator  Q(k)  is  stable  if  for  any  fixed  time  T  there 
exists  a  constant  Mr  such  that 

||Q"(fc)||  <  MT  (3.2) 

for  all  k  sufficiently  small  (say  k  <  ko)  and  nk  <  T. 

The  condition  (3.2)  ensures  that  for  all  initial  vectors  U°,  the  solution  Un  = 
Qn(k)U°  satisfies 

\\Un\\  <  Mt||£/°||  (3.3) 

for  nk  <  T. 

Iferc  |(  ■  |(  represents  some  norm  over  all  meshpoints  at  a  fixed  time.  For  example, 
the  discrete  1 2  norm  is  given  by  • 

\\un\\l  =  h  £  lie  I2 

m=-oo 

with  |  •  |  representing  the  usual  vector  two- norm. 

Up  until  Section  3.4,  whore  we  introduce  Sobolev  norms,  we  will  always  suppose  that 
the  norm  ||  ||  is  equivalent  to  the  £3  norm  in  the  sense  that  there  exist  constants  My  and 
M-i  such  that 

A/,||C/|ia  <  ||l/||  <  A/2III/H2 

Tor  all  U.  With  this  restriction,  Stability  Definition  3.1  is  independent  of  the  norm  used. 
If  Q(k)  is  stable  in  the  l?  norm  then  it  is  also  stable  in  any  equivalent  norm. 

The  following  equivalent  definition  of  stability  is  sometimes  easier  to  work  with 
since  it  only  requires  a  bound  on  ||Q(fc)J)  rather  than  a  uniform  bound  on  ||Qn(fc)||.  The 
dilficulty  in  applying  the  new  definition  is  that  such  a  bound,  when  it  holds,  will  often 
hold  only  in  a  very  special  norm  tailored  to  the  problem,  and  will  generally  not  hold  in 
equivalent  norms. 
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STABILITY  DEFINITION  3.1'.  The  operator  Q(k)  is  stable  if  there  exists  a  norm  ||  •  || 
and  a  constant  a  >  0  such  that 

IIQ(*)II  <  1  +  a*  (3.4) 

for  all  k  <  £<>• 

Clearly  (3.4)  implies  (3.2)  since 

\\Qnm  <  n<2(*)ir  <  (i + o*r  <  e°T 

if  nk  <  T  and  we  can  thus  take  Mr  =  eaT.  The  converse,  that  such  a  norm  exists  for 
any  stable  scheme,  is  proved  constructively  in  Chapter  4  of  Richtmyer  and  Morton[46] 
as  part  of  the  Kreiss  Matrix  Theorem. 

In  some  cases  bounds  of  the  form  (3.4)  can  be  obtained  directly.  This  method 
of  proving  stability  is  referred  to  as  the  energy  method  since  for  physical  systems  the 
required  norm  is  often  simply  the  energy  of  the  system.  Often,  however,  it  is  easier  to 
determine  stability  by  an  alternative  approach  known  as  the  von  Neumann  method.  We 
take  Un  to  be  a  single  Fourier  mode,  =  el(mhU  where  U  is  the  vector  of  Fourier 
coefficients  at  time  n,  and  insert  this  into  (3.1).  We  find  that  Un+X  is  again  a  single 
Fourier  mode  with  coefficients 

Un+i  =GU,k)Un 

for  some  matrix  C(£,  k),  called  the  amplification  matrix.  Stability  Definition  3.1  is 
equivalent  to  the  following  definition  based  on  this  amplification  matrix. 

STABILITY  Definition  3.2.  The  operator  Q(k)  is  stable  if  for  any  lixed  time  T 
there  exists  a  constant  Mr  such  that  powers  of  the  corresponding  amplification  matrix 
arc  uniformly  bounded  by  Mj ■, 

ll<?n(£,*)||  <  mt  (3-5) 

for  all  £,  k  <  k0  and  nk  <  T. 

Corresponding  to  Stability  Definition  3. 1'  we  have  the  following  definition  of  stability, 
which  is  again  equivalent. 

STABILITY  Definition  3.2'.  The  operator  Q(k)  is  stable  if  there  exists  a  norm  ||  •  || 
and  a  constant  a  >  0  such  that 

\\C({,k)\\<  l+ak  (3.6) 

for  all  (  and  k  <  fc0. 

Since  every  matrix  norm  is  bounded  below  by  the  spectral  radius,  we  find  from 
Stability  Definition  3.2'  that  a  necessary  condition  for  stability  is  the  so-called  von 
Neumann  condition : 

p(<7(t,*))  <  I  +  O(k).  (3.7) 

This  condition  is  frequently  sufficient  as  well.  Chapter  4  of  Itichtmyer  and  Morfon[46]  has 
a  thorough  discussion  of  sufficient  conditions.  Here  we  will  mention  ordy  a  lew  examples 
which  will  prove  particularly  useful. 

If  for  all  £  and  k,  C(£,k)  is  a  normal  matrix,  i.c.,  if  G  commutes  with  its  conjugate 
transpose,  then  ||(7(£,  fc)||y  =  p(G((,k)).  Ily  using  the  2- norm  in  Stability  Definition  3.2' 
we  sec  that  in  this  case  the  von  Neumann  condition  is  sufficient  for  stability. 

More  generally,  if  suffices  that  the  matrices  G((,k)  be  simultaneously  normalizable, 
as  defined  in  the  following  theorem  (see  Wichtmycr  Sr.  Morton(46]). 


THEOREM  3.1.  Suppose  there  exists  a  constant  matrix  S  such  that  SG(£,  k)S~l 
is  a  normal  matrix  Cor  all  £,  k  <  kg.  Then  the  von  Neumann  condition  is  sufficient  for 
stability. 

Proof.  Define  the  vector  norm  ||  •  ||s  by 

llvlls  =  WSyfo. 

This  vector  norm  is  equivalent  to  the  2-norm.  The  corresponding  matrix  norm  is 

llvlls  =  ||^5-,||2.  (3.8) 

In  this  norm  we  have 

(|G(e,fc)||s  =  ||5G(^*)5-1||2 
=  p(SGU,fc)S-‘) 

=  P(G((,*)) 

and  the  theorem  follows  by  using  the  norm  j|  ■  ||s  in  Stability  Definition  3.2'.  | 

An  important  application  of  this  theorem  provides  the  result  that  the  von  Neumann 
condition  is  sufficient  for  stability  if  the  G(£,  k)  are  simultaneously  diagonalizable  (since 
any  diagonal  matrix  is  normal).  Many  methods  for  the  problem  ut  =  Aux  have  the 
property  that  their  amplification  matrices  are  polynomials  in  the  matrix  A  and  hence 
arc  diagonalized  (for  all  £  and  A:)  by  the  eigenvector  matrix  of  A  (by  the  assumption 
of  hyperbolicity,  A  is  diagonalizable).  In  particular,  the  Lax-Wendroff  operator  and  the 
exact  solution  operator  have  this  property,  and  the  von  Neumann  condition  is  sufficient 
for  their  stability. 

3.2.  Stability  of  the  time-split  method. 

We  now  turn  to  the  stability  analysis  of  the  time-split  method  (1.23).  When  Q}(k/ 2)  = 
Qf(k),  as  is  true  for  the  splittings  (2.18),  for  example,  Cauchy  stability  of  the  Strang 
splitting  (1.21)  is  equivalent  to  stability  of  the  first  order  splitting 

t/n+l  =  Qf(k)Qa(k)Un.  (3.9) 

For  simplicity  we  restrict  our  attention  to  this  splitting,  and  set  Q(k)  =  Qf{k)Qa(k). 

Let  Gf(£,  k)  and  Ga(£,  k)  bo  the  amplification  matrices  corresponding  to  the  operators 
Qj{k)  and  ("^(A:)  ,  respectively.  Then  it  is  easy  to  verify  that  the  amplification  matrix 
G(£,k)  for  Q(k)  satisfies 

C(£,k)  =  Gf(£,k)G,(£,k). 

This  allows  us  to  calculate  the  amplification  matrix  for  the  time-split  method  relatively 
easily.  In  general  the  stability  of  Q/{k)  and  Qa(k)  separately  does  not  imply  that  f,>( fc)  is 
stable,  or  even  that  the  von  Neumann  necessary  condition  is  satisfied  for  G(£,  k),  since  the 
spectral  radius  is  not  submultiplicative  (i.c.,  the  inequality  p(G)  <  p(G f)p(Ga)  does  not 
hold).  It  is  easy  to  find  examples  for  which  Qf[k)  and  Qa[k)  arc  both  stable  operators 
but  (3.9)  is  unstable.  In  fact,  this  can  happen  even  when  Q/[k)  and  Qa(k)  arc  exact 
solution  operators  for  well-posed  hyperbolic  problems,  as  the  following  example  shows. 
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FIG.  3.1.  Spectral  radius  of  the  amplification  matrix  G(£,k)  of  Example  3.1 
for  p  =  5, 10,  as  a  funciton  of  ( k  between  —n  and  x. 


(Incidentally,  the  converse  can  also  occur,  i.e.,  the  product  may  be  stable  even  if  one  of 
the  operators  in  unstable  on  its  own.  See  Abarbancl  &  Gottlieb[l]  for  an  example  of 
such  a  scheme.) 

Example  3.1.  Let 


Then  the  problems  ut  =  A/ux  and  ut  =  Aaux  are  well-posed,  strictly  hyperbolic 
problems  for  any  value  of  the  parameter  p,  and  so  is  ut  =  (A f  +  A,)ux  if  p  >  —2. 
Let 

Qf(k)  =  exp(kA/dx),  Qa(k )  =  exp(kAadx). 

The  corresponding  amplification  matrices  are 

<?/(£,  k )  =  exp(j'<£/t f) 

_  e‘*f  pi  sin  k( 

.  0  e~xk*  . 


GM,k)  —  exp  (ik(Aa) 

cosk£  i  sin  k£ 
i  sin  k£  cos  k£  . 

We  have  p(Gf{£,k))  =  p(Ga(£,  k))  =  1  for  all  £  and  k.  On  the  other  hand,  the 
amplification  matrix  G((,k)  for  the  time-split  method  has  p{G(£,k))  =  1  for  all  (  and  k 
only  if  |/t|  <  2.  When  \p\  >  2,  the  method  (3.U)  is  unstable.  Figure  3.1  shows  graphs  of 
p(G(£,k))  for  p  =  5  and  10. 
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3.3.  Simultaneously  normalizable  splittings. 

As  wc  have  just  seen,  the  individual  stability  of  Qf(k)  and  Qa(k)  is  not  sufficient 
to  guarantee  the  stability  of  Q(k )  in  general.  However,  for  certain  special  cases  (which 
include  some  fairly  broad  and  important  classes  of  problems),  the  individual  stability 
is  sufficient  for  overall  stability.  Since  the  matrices  (?/(£,  k)  and  Ga((,k)  are  generally 
much  easier  to-work  vv.Ui  than  their  product  C((,k),  it  is  useful  to  identify  such  classes 
of  problems.  For  these  problems  stability  is  relatively  easy  to  determine. 


We  first  note  that  if  there  exists  a  norm  ||  ■ 

||  and  a  constant  a  such  that 

||G/(£,*)||  <  1+afc 

V  £,  k  <  k0 

(3.10a) 

and 

\\Ga(t,k)\\  <  l  +  <xk 

V  £,  k  <  k0. 

(3.10b) 

\m,k)\\  <  ||C/(*,*)||||C.tf,*)ll 

<  1  +  2  ak  +  a2k 2 

<  t  +  a0k  V  £,  k  <  k0 
where  ao  =  2a  +  a2fco,  so  Q(k)  is  stable. 

Of  course  if  Qf{k)  is  stable  then  by  Stability  Definition  3.2'  there  exists  a  norm  such 
that  (3.10a)  holds.  Similarly,  if  Qa(k)  is  stable  then  (3.10b)  also  holds  in  some  (possibly 
different)  norm.  Only  in  certain  special  cases  can  wc  easily  show  the  existence  of  a  single 
norm  in  which  both  (3.10a)  and  (3.10b)  hold. 

As  one  such  case,  suppose  that  all  of  the  matrices  Gf(£,k)  and  (?,(£,  k)  are  simul¬ 
taneously  normalizable  by  a  single  matrix  S.  Then  the  individual  stability  of  Q/(k)  and 
Qa(k)  guarantees  the  stability  of  Q(k),  since  then  (3.10a)  and  (3.10b)  both  hold  in  the 
5-norm  defined  in  (3.8). 

Some  operators,  such  as  LW  and  exact  solution  operators,  have  the  property  that 
if  the  coefficient  matrix  is  normal  then  the  corresponding  amplification  matrix  will  also 
be  normal,  for  all  £  and  k.  Restricting  our  attention  to  such  schemes,  wc  find  that  it 
then  suffices  for  the  stability  of  Q(k)  that  the  two  matrices  Af  and  A„  be  simultaneously 
normalizable  and  that  Qf(k)  and  Qa(k)  be  individually  stable. 

This  result  is  quite  useful,  since  in  many  practical  problems  the  matrices  A/  and 
Aa  are  simultaneously  normalizable.  This  class  includes,  for  example,  scalars,  symmetric 
matrices,  and  commuting  matrices  (which  are  simultaneously  diagonalizablc). 

These  results  can  easily  be  extended  to  splittings  involving  more  than  two  terms. 
Since  this  is  frequently  useful,  wc  summarize  the  above  results  and  their  proofs  in  a  more 
general  setting. 

TllKOItHM  3.2.  Let  A\,  A-2, . . .,  Am  be  constant  matrices.  Approximate  each  solution 
operator  ex\)(kA}Ox)  by  some  operator  Qj(k)  with  amplification  matrix  6?y(f,  k).  Suppose 
there  exists  a  single  norm  j|  •  |j  and  a  constant  a  such  that 

l|Gjtt.*)ll  <  1+a*  V{,  k  <  k0,  j=  l,2,...,m.  (3.11) 

Then  the  scheme 

Un+l  =Qi(kl)Q,(k'i)..-Qm{km)Un  (3.12) 
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is  stable. 

Proor.  Let  Q(k)  =  <?,(*)•••£«(*)  and  let  G[(,k)  =  Ct((,k)  -  Cm(£,  k)  be  the 
corresponding  amplification  matrix.  Then 

iiG(e,fc)ii<iiG1u,fc)ii-.-iiGm(e,fc)ii 

^  1  +  oo^ 

for  ao  =  ma  +  (”)asJfco  +  •  •  •  +  and  hence  Q(fc)  is  stable.  | 

Theorem  3.3.  With  the  A}and  Gj((,  k)  as  in  Theorem  3.2,  suppose  there  exists 
some  nonsingular  matrix  S  such  that  SGj(£,  k)S~l  is  a  normal  matrix  for  all  j,  and 
k.  Suppose  furthermore  that  each  satisfies  the  von  Neumann  condition, 

p(Gj((,  k))  <  1  +  ak  V  £,  *  <  *o,  j  =  1, 2, . . .,  m 

for  some  constant  a.  Then  Q(k)  is  stable. 

Proof.  Using  the  5- norm  defined  in  (3.8)  and  the  fact  that  SGj(£,  k)S~l  is  normal, 
wc  h^vc 

l|G>(£,A)iis  =  HSCitt,  fcjs-'ii* 

=  p(SGj(i,k)S-') 

=  p(G](i,k)) 

<  l  +  ak 

and  stability  follows  by  Theorem  3.2.  | 


3.4.  Block  triangular  systems. 

A  similar  stability  result  can  he  obtained  for  the  standard  block  triangular  system 


u 

II 

id 

u 

V 

i  l  0  Am 

V 

with  the  splitting  (1.32).  The  solution  v  does  not  depend  on  u.  In  solving  for  u,  the 
computed  vx  enters  essentially  as  a  forcing  function.  Because  of  this  we  obtain  only 
a  weak  stability  result,  in  wheih  the  norm  of  ||(/n||  is  bounded  in  terms  of  a  discrete 
Sobolev  norm  of  the  initial  conditions.  The  Sobolev  .iorm  |||  •  ))|  has  the  form 

mi=iif/n+iiD+f/||. 

With  the  splitting  (1.32),  the  schemes  Qa(k)  and  Q/(k)  will  be  of  the  form 


Q.(k)  = 


I 

.0 


Qn{k)' 

CM*)]’ 


<?/(*)  = 


Qu(k)  O' 
0  / 


(3.13) 


Suppose  that  Qu(fc)  and  Qn(k)  arc  stable  schemes.  Then,  in  particular,  there  exists  a 
norm  j|  •  ||  and  a  constant  a  >  0  such  that 


||<?i,(fc)||<  1+afc  V*<*0.  (3-14) 
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All  or  the  following  estimates  will  be  in  this  norm.  We  also  suppose  that 


||<3«(*)V||  <  kM\\D+V\\  VV,k<k0  (3.15) 

for  some  constant  M.  For  example,  if  Q,(k)  =  LW(At,  k),  we  have 

#  s 

Qia(k)  —  kAi^Do  +  jfc2  Aia  A22D+D— 

—  ^kAn{D+  +  i)_)  +  —  Z)_) 

since  D+D _  =  (D+  —  D-)/h.  Since  ||Z5_ Vjl  =  (ID+KH,  we  have 

lieta(*)VH  <  k(HAi2ll  +  X||A,2A22||)||Z>+K||. 

For  fixed  X  this  is  of  the  form  (3.15). 

With  these  assumptions  we  then  have  the  following  theorem. 

THEOREM  3.4.  Suppose  Q/(k)  and  Qa(k)  are  stable  schemes  as  above.  Then  the 
split  scheme  Qa(k)Qj(k)  is  weakly  stable: 

||f/»||  <  KT(\\U°\\  +  ||D+V°||)  (3.16a) 

||V"||  <  kT\\V°\\  (3.16b) 

for  nk  <  T.  Here  Kt  and  Kt  are  constants  depending  only  on  the  fixed  time  T. 

Proof.  When  the  full  scheme  f/"+l  =  Q a[k)Q f[k)Un  is  written  out  we  obtain 

t/n+1  =  Qn[k)Un  +  Qii(k)Qla(k)Vn  (3.17a) 

Vrn+l  =  Q22{k)Vn.  (3.17b) 

The  bound  (3.16b)  follows  immediately  from  (3.17b)  and  the  stability  of  Q22(k).  Moreover, 
by  linearity,  an  identical  bound  holds  for  the  linear  combination  of  solutions  D+Vn,  i.c., 

\\p+vn\\  <  kT\\D+v°\\. 

Using  this  together  with  (3.15)  in  (3.17a)  gives 

l|f/nM||  <  HOn(*)ll(ll^/n|l  +  *AfA('r|(£>+K0||). 

When  iterated  n  times  this  gives 

\wn\\  <  iiQ..(*)iriit/°ii  +  +  iiQ,,(fc)ir-a 

V  (3.18) 

+  -"  +  ll«u(fc)||  +  lJ||/>+V°||. 

By  (3.14),  i|Qit(*)ir  <  (1  +  «*)"  <  eaT  ir  nk  <  T.  Using  this  in  (3.18)  gives 

\\Un\\  <  eaT(\\U0\\  +  TMkT\\D+V°\\) 
for  nk  <  T,  which  is  of  the  desired  form  (3.16a).  | 
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3.5.  The  shallow  water  equations. 

We  will  now  investigate  the  stability  of  the  splitting  (1.38)  for  the  shallow  water 
equations.  Since  this  is  a  quasi  linear  system  of  equations,  a  complete  stability  analysis  is 
diflicuit  to  perform,  even  for  unsplit  methods.  We  will  perform  only  a  linearized  stability 
analysis  for  the  corresponding  frozen  coefficient  problem  with 


Af  =  - 


®  0o/2 

A  —  _ 

U0  T0 

0  o/2  0 

i  **9  — 

To  Uq 

(3.19) 


Here  the  constant  Uq  is  a  representative  value  of  u  while  To  is  a  representative  value  of 
(0  —  ^0)/2.  One  hopes  that  if  a  method  is  stable  on  the  frozen  coefficient  problem  for  all 
values  of  Uq  and  T0  in  the  appropriate  range,  then  the  method  will  also  be  stable  on  the 
nonlinear  problem.  It  is  well  known  that  this  is  not  necessarily  so;  nonlinear  instabilities 
may  arise.  Nonetheless,  the  linearized  stability  analysis  is  valuable  because  an  instability 
for  the  frozen  coefficient  problem  will  almost  certainly  lead  to  instability  of  the  nonlinear 
problem,  and  thus  we  at  least  obtain  upper  bounds  on  the  stability  limit.  Moreover,  for 
the  shallow  water  equations  computations  indicate  that  the  nonlinear  scheme  is  usually 
stable  when  the  frozen  coefficient  problems  are. 

Stability  of  the  scheme  (2.42)  applied  to  (3.19)  is  easy  to  determine  using  Theorem 
3.3.  The  matrices  A/  and  A„  arc  both  symmetric  and  so  Q{k )  is  stable  provided  Qf(k) 
and  Q,(k)  are  both  stable.  Since  Qf(k)  is  the  exact  solution  operator,  it  is  always  stable, 
and  so  stability  is  determined  entirely  by  Q3(k).  The  eigenvalues  of  A,  arc  U0  ±  To  and 
so  Qa(k)  —  LW(A„,  k)  is  stable  if  (Uq  ±  To )k/h  <  1. 

Suppose  that  (1.37)  holds,  i.c.,  |u|  <  f0o  and  |0  —  0o|/2  <  e0o  for  all  x  and  t  for 
the  solutions  of  interest.  Then  all  of  the  relevant  frozen  coefficient  problems  are  stable 
provided 


h  2«f>o- 


(3.20) 


Note  that  for  the  unsplit  method  LW(A,  k),  the  stability  limit  is  roughly 


The  split  scheme  is  thus  stable  for  much  larger  values  of  k.  Recall,  however,  from  Section 
2.9  that  for  the  split  method  an  accurate  solution  is  obtained  most  efficiently  using 
k/h  sa  l/0o-  Such  mesh  ratios  are  well  within  the  stability  limit  (3.20)  and  long-time 
calculations  on  the  full  nonlinear  system  have  revealed  no  instabilities. 


54 


4.  Boundary  conditions  for  the  intermediate  solutions 


4.1.  Introduction  and  a  simple  example. 

So  Tar  wc  have  considered  the  time-split  method  applied  only  to  the  Cauchy  problem 
on  the  unbounded  spatial  domain  or  to  problems  with  periodic  boundary  conditions. 
In  practice  we  must  be  able  to  deal  with  more  general  boundary  conditions.  The 
implementation  of  finite  difference  schemes  frequently  requires  more  boundary  data  than 
arc  supplied  with  the  differential  equation.  In  particular,  when  using  a  time-split  method, 
special  boundary  data  must  be  generated  for  the  intermediate  solutions. 

For  the  most  part  wc  will  restrict  our  attention  to  the  time-split  method  (2.18a,b) 
for  solving  perturbed  hyperbolic  problems,  although  the  same  techniques  can  be  applied 
to  a  wide  variety  of  other  problems  and  splittings.  Some  examples  of  other  applications 
are  given  in  Chapter  5. 

We  begin  our  discussion  with  a  simple  example  which  illustrates  the  problems 
encountered  and  the  general  methodology  used  to  determine  the  correct  boundary  data. 

A  constant  coefficient  scalar  problem.  Consider  the  equation 


«t  —  -(1  +  e)u* 

(4.1) 

on  the  strip  0  <  x  <  1,  t  >  0,  with  initial  conditions 

u(x,0)  =  /(x),  0  <  x  <  1, 

(4.2) 

and  boundary  conditions 

«(0,  <)  =  g(<),  t  >  0. 

(4.3) 

For  c  >  — 1,  this  is  a  well-posed  problem  as  it  stands.  Boundary  data  is  prescribed  only 
at  the  inflow  boundary  x  =  0.  Values  at  the  outflow  boundary  x  =  1  are  determined  as 
part  of  the  solution. 

The  exact  solution  to  this  problem  is  a  wave  moving  to  the  right,  unaltered,  with 
speed  l  +  c: 

u(x,<)  =  J(x-  (l  +  <)*),  0  <  x  <  1,  t  >  0 

where  for  £  <  0  wc  define 


/(fl  =  »(-*/(  1+0),  £<0. 


We  will  first  consider  the  unsplit  Lax-WcndrolT  method.  If  the  mesh  spacing  in 
the  x-direclion  is  h  —  I / N  for  some  integer  N,  then  the  grid  points  of  interest  arc 
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xo,xi,...,xn.  The  Lax-Wendroff  method  is 

=  +  U4) 

+  ixa(l  +  c)2(t/~ +1  -  2 IC  +  «C-i).  m  =  1,2, . . iV  —  1.  1  •  1 

m  *» 

This  scheme  cannot  be  applied  (or  m  =  0  or  m  =  N  and  so  f/ J+1  and  t/Jy+1  must 
be  determined  in  some  different  manner.  At  the  left  boundary  we  simply  use  the  given 
boundary  data  (4.3), 

t/o+‘  =  9(tn+x).  (4.5) 

At  the  outflow  boundary  we  must  cither  extrapolate  from  the  interior,  e.g., 

UnN+1  =  2U"n+_\  -  UtfJ*,  (4.6) 


or  use  a  one-sided  difference  scheme,  e.g., 

W  =unN-  X{1  +  ()(irN  -  U%_x).  (4.7) 

Both  (4.6)  and  (4.7)  have  local  truncation  errors  which  are  0(k2).  This  is  sufficient  to 
retain  the  0(k 2)  global  accuracy  of  the  Lax-Wendroff  method.  In  general  the  overall 
accuracy  of  a  method  is  not  degraded  by  errors  in  the  boundary  values  provided  the 
local  error  at  the  boundary  is  no  larger  than  the  global  error  for  the  interior  scheme  (see 
Gustafsson[28|).  In  addition,  of  course,  the  total  method  (including  boundary  schemes) 
must  be  stable.  Stability  is  more  difficult  to  determine  for  initial  boundary  value  problems 
than  for  Cauchy  problems  and  is  discussed  in  Section  4.6.  For  this  simple  problem  both 
(4.6)  and  (4.7)  yield  stable  methods. 

Now  consider  a  time-split  method  applied  to  the  same  problem  (4.1)  with 

Af  =  -1,  A,  =  c. 

We  now  assume  that  eCl.  Since  the  operators  commute,  there  is  no  need  to  use  the 
Strang  splitting  and  so  wc  need  introduce  only  one  intermediate  solution.  Taking  k  =  ph 
for  some  integer  p  >  I  and  using  the  exact  solution  operator  on  the  fast  part  together 
with  1,W(A„  k),  the  split  method  is 

V*m  =  um-P<  m  =  p,p+  +  1, 

unm+l  =  u‘m-  "  tC-i)  +  ipVttC-M  "  2 1C  +  lC_,), 

m  —  1,2 

Notice  that  we  use  (4.8a)  to  define  f/Jv+i  even  though  it  is  not  within  the  domain  of 
interest.  Nonetheless,  it  can  be  used  in  computing  UJf''  (which  is  of  interest)  in  the 
Lax-Wendroff  step  (4.8b).  Because  or  this  we  do  not  need  any  special  procedure  to  specify 
I/n+i.  This  is  one  advantage  of  using  time-split  methods  for  such  perturbed  problems. 
Since  they  arc  essentially  skewed  (one-sided)  Lax-Wendroff  methods  which  follow  the 
characteristics  of  the  problem,  artificial  boundary  values  are  often  not  needed  at  outflow 
boundaries. 


(4.8a) 

(4.8b) 
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Instead,  we  need  to  specify  additional  values  at  x  =  0.  Wc  still  use  (4.5)  for  U q+1, 
but  wc  must  also  specify  U0,Ul,...,Up_i.  Alternatively,  wc  can  leave  these  values 
unspecified  and  apply  (4.8b)  only  for  m  =  p  +  We  must  then  determine 

f/"+1,...,t/p+1  by  some  alternative  procedure. 

Since  there  is  no  splitting  error  for  this  problem,  the  results  of  Section  2.5  indicate 
that  for  optimal  efficiency  we  should  take  p  1/e.  However,  for  simplicity  we  first 
consider  the  case  p  =  1.  Then  we  only  need  to  specify  Uq  or  l/"+1. 

Three  possibilities  for  specifying  U"+1  are  immediately  apparent.  The  first  is  to 
interpolate  between  the  known  values  U j}+1  and  U J+1, 

f/r+1  =  &US+1  +  t/J+1).  (4.9) 

This  is  0(k2)  accurate.  However,  when  e  is  small  this  choice  causes  a  severe  loss  of 
accuracy  in  (4.8)  and  completely  negates  the  increase  in  efficiency  obtainable  through 
the  use  of  the  time-split  method.  The  reason  is  that  the  local  truncation  error  for  the 
method  (4.8)  is  0(ek3)  giving  0(ek 2)  global  errors.  It  is  this  factor  of  t  which  makes  the 
time-split  method  advantageous  over  the  unsplit  method  (4.4).  Tly  using  (4.9)  wc  lose 
this  advantage. 

Figures  4.1a,b  show  the  errors  at  time  t  =  0.4  using  this  time-split  method  with  the 
boundary  conditions  (4.9)  when  e  =  0.1.  Signals  propagate  with  velocity  1  +c  =  1.1  and 
so  errors  from  the  improper  specification  of  U*+l  have  propagated  in  to  approximately 
x  =  0.44  at  this  time.  To  the  right  of  this  point  all  errors  are  due  solely  to  the  interior 
scheme.  It  is  this  accuracy  which  wc  would  like  to  match  at  the  boundary.  Clearly 
the  boundary  approximation  (4.9)  is  causing  a  loss  of  accuracy.  When  t  is  smaller,  as 
in  Figures  4.1c,d  where  t  =  0.001,  this  disparity  in  the  size  of  the  errors  is  even  more 
apparent. 

In  order  to  maintain  th*'  advantage  of  the  time-split  method,  wc  must  use  a  more 
accurate  boundary  scheme,  one  with  local  error  0(ek2).  One  possibility  is  to  use  higher 
order  interpolation.  Using  quadratic  interpolation  on  the  points  (/q+1,  U£+l,  6/3"*" 1 
would  give  0(k3)  errors.  For  k  sufficiently  small  ( k  <  e),  this  provides  sufficiently  ac¬ 
curate  data.  However,  th  '  use  of  higher  order  interpolation  can  cause  stability  problems. 
Moreover,  when  p  >  1  there  will  be  several  values  l/"+l, . . . ,  t/”+1  to  be  determined  and 
interpolation  is  unsatisfactory. 

The  second  obvious  choice  for  f/"  +  1  is  to  simply  use  Lax-Wendroff  on  the  unsplit 
problem, 

(/?+’  =  LW(-(1  +  0.  (4.10) 

This  also  has  0(fc3)  local  error  and  provides  sufficiently  accurate  data  for  small  k.  Again, 
however,  stability  may  be  a  problem  and  Tor  p  >  1  the  scheme  is  certainly  unstable. 

The  final  approach  to  specifying  '  is  based  on  Taylor  series  expansions  in  from 
the  boundary.  This  is  the  best  approach  and,  for  this  simple  problem,  gives  the  correct 
value  of  f/"+l  exactly.  We  want  f/"+1  to  approximate  u(h,tnj  1).  Wc  can  expand  this 
in  a  Taylor  scries  about  «(0,  tn+ 1): 

u{h,tn+ 1)  =  u{0,tn+i)  +  hux(Q,  fn+i)  +  j/»2Uix(0,tn+i)  +  •••.  (4.11) 

Approximating  this  directly  by  differencing  the  known  values  f/y  +  ’  would  give  us  the 
interpolation  scheme  rejected  above.  However,  using  the  differential  equation  (4.1)  we 
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(a)  k  =  1/25,  e  =  10~l 


(b)  k  =  1/100,  £  =  IQ"1 
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(c)  k  =  1/25,  £  =  10~* 
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FlG.  4.1.  Errors  in  the  computed  solution  ot  (4.1)  using  the  split  scheme 
(4.8)  with  p  =  1  and  the  interpolator;/  boundary  condition  (4.9).  The  errors 
arc  shown  on  a  logarithmic  scale  for  various  values  of  k  and  f  .  Note  that 
the  interior  error  is  0(<k 2)  while  the  boundary  error  is  0(ka). 
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ii  uu  mat 


3  >  0, 


(4.12) 


d{u 


Wc  can  thus  rewrite  (4.11)  as 


% 


“(^>  <n  +  l)  =  **(0»  tn+l)  ~  T+7u«(0»  *»  +  l)  +  Kl+«)  u«(®>  *»+»)  +  -  ‘ • 


(4.13) 


The  desired  data  is  now  expressed  in  terms  of  1-derivatives  of  u  along  the  boundary,  i.e., 
derivatives  of  the  known  function  g(t)  from  (4.3).  For  this  simple  problem  (4.13)  can  in 
fact  be  evaluated  in  closed  form,  giving  the  desired  value  of  U”+l  exactly: 

U?+' =  g(tn+t  -  h/( l+£».  (4.14) 


The  calculations  shown  in  Figure  4.1  have  been  repeated  using  (4.14)  instead  of  (4.9). 
The  results  are  shown  in  Figure  4.2.  Since  for  this  problem  the  boundary  data  (4.14)  is 
exact,  the  errors  are  actually  smaller  near  the  boundary  than  in  the  interior. 

This  same  approach  can  be  used  in  a  wide  variety  of  problems  to  determine  boundary 
data  for  points  near  the  boundary.  In  general  it  will  not  be  possible  to  obtain  the  exact 
data  in  closed  form  as  in  (4.14),  but  a  series  solution  can  be  developed  and  evaluated 
to  arbitrary  accuracy.  Goldberg  &  Tadrnor(21](22]  explain  how  to  do  this  for  general 
inflow-outflow  boundaries.  This  will  also  be  discussed  in  Section  4.4. 

It  seems  that  we  have  completely  avoided  the  need  to  specify  boundary  values  for 
the  intermediate  solution  U' .  For  this  simple  problem  that  is  true.  However,  for  many 
problems  it  is  not  possible  to  avoid  specifying  intermediate  boundary  values.  This  is 
particularly  true  when  implicit  methods  are  used  in  the  splitting.  In  other  situations  it  is 
simply  more  convenient  computationally  to  specify  boundary  values  for  the  intermediate 
solution  than  to  leave  these  points  unspecified. 

The  remainder  of  this  chapter  is  devoted  to  showing  how,  Tor  many  problems,  the 
same  approach  used  above  to  compute  f/”+l  can  be  extended  to  compute  arbitrarily 
accurate  intermediate  boundary  data. 

Computing  U*0.  We  now  return  to  our  original  plan  to  specify  U*0  Tor  the  scalar 
problem  (4.1).  We  require  data  at  the  point.  *o>  which  is  on  the  inflow  boundary.  At 
this  boundary  the  data  (4.3)  has  been  supplied,  but  is  not  usable  directly  since  U  is 
obtained  not  by  solving  the  original  equation  but  rather  by  solving  the  subproblcm 
«,*  =  —  tt*.  This  is  the  fundamental  step  in  correctly  computing  intermediate  boundary 
data:  introduce  a  new  fund  ton  u  which  solves  the  differential  equation  actually  being 
approximated  in  the  relevant  step  of  the  splitting.  The  desired  boundary  data  can  then 
be  expanded  as  a  Taylor  scries  in  this  function.  In  many  cases  this  can  be  reexpressed  as 
a  series  in  the  original  variable  and  evaluated  in  terms  or  g(l)  as  before.  We  will  sec  that 
for  many  problems  it  is  possible  to  generate  stable  0(<kz)  boundary  data  quite  easily. 
For  I  he  problem  (4.1)  we  can  in  fact  generate  boundary  data  which  is  exactly  correct, 
just  as  wc  did  for  U”+l . 

Consider  a  single  step  (4.8a)  of  the  lime-split  method  starting  at  time  tn  and  suppose 
that  f/JJ,  =  u(j'm,  /,,).  Since  in  this  problem  wo  have  used  the  exact  solution  operator 
exp(k/\ f<)z)  in  (4.8a),  U  is  then  the  exact  solution  at  time  tn  M  to  the  subproblcm 


«t  —  -u*  *>(>,<>  tn, 


(1.15) 
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k  =  1/25,  t  =  10 


—  in-* 


*  =  1/25,  £  =  10-* 
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Fig.  4.2.  Errors  in  the  computed  solution  or  (4.1)  using  the  split  scheme 
(4.8)  with  p  =  1  mid  the  correct  boundary  condition  (4.14).  The  errors  arc 
shown  on  a  logarithmic  scale  for  various  values  of  k  and  (. 


with  initial  conditions 


«*{M«)  =  *(M«)» 


x  >  0 


(4.16) 


at  time  tn.  The  idea  is  to  use  the  differential  equations  (4.1)  and  (4.15)  to  transform  the 
given  boundary  conditions  (4.3)  for  u  into  boundary  conditions  for  u*.  We  wish  to  find 
an  appropriate  value  for  UQ,  which  should  be  an  approximation  to  «‘(0,  tn+j).  This  we 
can  expand  in  a  Taylor  scries.  Using  (4.15),  we  find  that 


tt*(0,  tn  +  k)  —  u*(0,  tn)  +  fcu,*(0,  tn)  +  |*2ut*JO,  tn)  +  • 
=  u*(0,  tn)  -  ku*x( 0,  tn)  +  £*2tt*x(0,  tn)  + 


(4.17) 


Since  the  initial  conditions  (4.16)  hold  for  all  x,  that  relation  can  be  differentiated  with 
respect  to  x,  giving  u‘(x,tn)  =  ux(x,tn)  and  similarly  for  higher  derivatives.  So  (4.17) 
becomes 

tt*(0,f„  +  *)  =  u(0,tn)  ~  kvx(0,  <„)  +  |*2uXI(0,fn)  +  •••.  (4.18) 

We  can  now  use  the  original  equation  (4.1)  governing  u  to  rewrite  this  in  terms  of  t- 
derivatives  of  u.  Using  (4.12),  (4.18)  becomes 


U*(0,  tn  +  k)  =  u(0,  tn)  +  «,(0,  tn)  +  3(j^)2«tt(0,  tn)  + 

=  g(tn  +k/(\  +c)). 


(4.19) 


This  is  the  desired  boundary  data  f7y,  expressed  in  terms  of  the  given  boundary  data 
(4.3). 

For  such  a  simple  example  it  is  easy  to  verify  that  this  is  the  correct  boundary  value. 
According  to  the  scheme  (4.8a)  we  would  really  like 

Uo  =  U'. 1,  =u(-h,tn)- 


(Recall  that  p  =  t.)  Of  course  u  is  not  really  defined  for  x  <  0,  but  using  the  differential 
equation  (4.1)  it  can  easily  be  extended  backwards  in  time  from  the  boundary.  Since 
(4.1)  has  characteristics  with  slope  1/(1  +  t),  wc  find  that 

u(-~h,  tn)  =  u(0,  tn  +  h/(\  +  c))  =  g(tn  +  k/(l  +  r )) 
exactly  as  in  (4.19). 

Uecausc  the  characteristics  for  the  problems  (1.1)  and  (4.15)  have  different  slopes, 
we  see  that  the  value  u ( — /*,  £n)  is  equal  to  both  «((),  tn  +  k/0  +  r))  and  u*(0,  t„  +  k)  and 
therefore  they  are  equal  to  each  other.  This  is  illustrated  in  Figure  4.3. 

Whim  p  >  l  we  can  compute  l/;  for  0  <  j  <  p  in  a  similar  manner.  Using  the 
fact  that  we  know  the  exact  solution  operator  for  the  subproblcm  (4. 15),  wc  can  project 
these  values  back  to  the  boundary  along  the  characteristics, 

f/*  =  u'(jh,tn  +  k) 

=  «*((),  t„  +  k-jh). 
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i  =  0 


1  his  boundary  value  can  bo  computed  as  before,  giving  the  general  expression 

Vj  =  9(t»  +  (k- jh)/(  l +()),  j  =  0, 1.  (4.20) 

Compulations  confirm  that  the  use  or  this  boundary  value  Tor  U"Q  in  the  split  method 
(4.8)  gives  excellent  results  that  are  virtually  identical  to  those  seen  in  Figure  4.2. 

In  general  when  using  this  approach  to  specify  boundary  conditions  for  the  inter¬ 
mediate  solutions  it  will  not  be  possible  to  generate  exact  boundary  data  as  we  did  here. 
It  often  will  be  possible,  however,  to  develop  a  scries  solution,  as  in  the  first  line  of 
(4.10),  which  can  be  used  to  generate  arbitrarily  accurate  boundary  data.  In  the  next 
lew  sections  we  demonstrate  how  this  can  be  done  Tor  systems  of  increasing  complexity, 
culminating  in  Section  4.5  with  the  development  of  boundary  conditions  for  the  shallow 
walcr  equations,  a  quasilinear  system  of  equations  with  inllow-outllow  boundaries. 

4.2.  Constant  coefficient  systems — inflow  boundaries. 

As  the  next  step  in  this  direction,  consider  a  constant  coefficient  system  of  equations 

«t  =  Aux  =  (A/  4-  A,)uxi  x  >  0,  t  >  0, 

“(*,0)  =  /(*),  (4.21) 

«(0,  t)  =  g(t),  t  >  0, 
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on  the  quarter  plane  x, t  >  0.  We  assume  that  the  boundary  x  =  0  is  a  pure  inflow 
boundary,  i.c.,  that  A  has  strictly  negative  eigenvalues.  We  also  assume  that  A/  has 
nonpositive  eigenvalues.  In  general  A/  and  A ,  do  not  commute,  so  we  will  have  to  use  a 
Strang-type  splitting.  There  will  be  at  least  two  intermediate  solutions,  say 

U-^cxpUkAfdjU" 

U **  «  e\p(kAMds) cxp(^kAfdx)Un . 

Of  course  there  may  be  many  more  if  exp(^kAfdx)  is  itself  approximated  by  several 
steps  of  Lax-Wendroff,  but  they  can  be  handled  similarly.  The  general  principle  should 
be  clear  from  considering  (4.22). 

Again  introduce  the  function  u(x,  t)  which  satisfies  the  first  subproblem  of  interest, 

u’  =  Afux,  x  >  0,  t  >  tn, 

u*(x,tn)  =  u(x,tn),  x  >  0. 

We  then  want 

U0  —  u  (0,  <„+ 1/2) 

=  «*(0,  t„)+  <„)  +  |fc2ut*t(0,fn)  +  ••• 

=  «*(0,  <„)  +  ;|A:A/«‘(0,  tn)  +  0,  <„)  +  ••• 

=  u(0,  tn)  +  %kAfuz[ 0,  tn)  +  \k2A)uxx{ 0,  <„)  +  ••• 


(4.23a) 

(4.23b) 


(4.24) 


where  we  have  used  (4.23a)  to  replace  ^-derivatives  of  u  by  z-derivatives.  I'll  esc  were 
then  replaced  with  z-dc.viva.Uvcs  of  it  using  the  initial  conditions  (4.23b).  We  next  use  the 
original  equation  (4.21)  to  replace  z-derivatives  of  u  by  ^-derivatives,  which  are  equivalent 
to  derivatives  of  the  boundary  data, 

U'a  =  n(0,tn)+'tkAfA-lut(0,tn)+yk2A'}A-2uu(0,tn)  +  --- 

=  g(tn)  +  ±kAfA-' +  tk^A'jA  -2<f{tn)  +  ■■:  (  ‘  °} 

We  have  assumed  that  A  has  strictly  negative  eigenvalues  and  thus  is  invertible.  In 
general  f/„  must  now  be  approximated  by  the  first  lew  terms  of  (4.25).  Keeping  the  first 
three  terms  gives  ()[k  ’)  accurate  boundary  data.  As  usual,  this  is  sufficiently  accurate  if 
k  is  small.  However,  it  is  worth  pointing  out  that  we  can  frequently  achieve  the  (){<k2) 
accuracy  we  desire  more  easily.  Suppose  ||/t",|(  =  0(1).  Then  since  A  =  Af  +  ()(<), 


A3fA  3  =  I  +  ()(< )  for  j  —  1,2,.... 

We  can  then  retain  0(ik2)  accuracy  simply  by  taking 

Vo  =  <j(tn  >  I/*)  +  WA,A-'  -  l)‘/(tn).  (4.20) 

We  may  still  wish  to  use  additional  terms  of  the  expansion  in  order  to  ensure  that  the 
error  from  the  boundary  conditions  does  not  dominate  the  interior  error.  The  boundary 
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conditions  (4.26)  arc  the  correct  order  of  accuracy  but  the  error  constant  may  be  larger 
than  that  of  the  interior  scheme.  We  obtain  0(ck3)  accurate  boundary  data  by  using 

Ul  =  g(<„+1/2)  +  WAfA~'  -  l)g’(tn)  +  -  I)g"(tn).  (4.27) 

The  additional  Work  incurred  by  using  three  terms  of  the  expansion  rather  than  two  at 
the  boundary  is  negligible  compared  to  the  work  being  done  in  the  interior. 

Now  to  find  boundary  values  for  U *  .  The  easiest  way  to  proceed  is  to  note  that 

U**  =  exp(-  ^kA/dx)Un+i 

which  prompts  us  to  define  u**(i,t)  as  the  continuous  solution  to 


wt**(a;,  t)  =  Afu“(x,  t)  x  >  0,  t  <  t„+l 
u‘*(x,tn+l)  =  u(x,t„+l)  x  >  0. 

We  now  solve  this  backwards  in  time  for 

U0  =«  (<Mn+ 1/2)- 

I’roceeding  as  in  (4.24)  and  (4.25)  we  obtain 

^0*  =  9{tn+ 1)~  hkAfA~lg'{tn+\)  +  ^k2A2fA~2g''{tn+i)  + 
+ 1/2)  -  {k{AfA~x  -  l)g'{tn+l). 


(4.28) 


Example  4.1  Consider 


u 

[-l 

Cl 

V 

1; 

e 

,f2 

—2 

y 

0  <  x  <  1,  t  >  0, 


w(:c,  0)  =  f(x),  0  <  x  <  1, 

il(0,  l)  —  g(t),  t>  0, 


where  u  —  (•«,  v)T .  For  the  splitting  we  take 

Af  = 


-1  0 
0  -2 


= 


0  c,l 

Ua  0 


Using  (2.22)  the  splitting  error  is  computed  to  be 


Ht(fc)  =  - 


L  4°2  fi « 2 


k9i. 


If  we  use  the  time-split  method  (2.l8a,b)  then,  according  to  (2.24),  the  optimal  stepsizo 
ratio  is 


U  4-  f ; 
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where  c  =  max|fy|.  For  Jfe  =  2h  and  h=  l/N,  (2.18a, b)  becomes 

=  rn  =  1,2 . N 

^  =  ^-2.  m=  2,3,  ...,1V 

OZ  =LW{A„k)0*m,  m  =  1,2 . IV -1 

t>o+l  =  »(*-+0 

tC+‘ =*C- 1.  m  =  l,2,...,lV 

Vl+l  =  V”_a,  m  =  2,3,...,fV. 

Notice  that  no  boundary  conditions  need  to  be  specified  at  the  outflow  boundary  z  =  1. 
On  the  inflow  ride  we  still  need  to  specify  UQ,  V],  UQ  ,  and  V7+l.  For  this  problem, 

.a._a _ 1  4  +  C|fa  3ei 

^  (2  — £i£j)2  12£2  4  +  4£i£j 

=  I+0(e). 

and  we  can  retain  0(cA;3)  accuracy  aking 

U0  =  ff(*«+i/a)  +  5 k(A/A~l  -  I)rf{tn) 

Similarly  we  use 

U'Z  =  9(tn+i„)  -  \k(A,A~'  -  I)g\tn+l). 

Wc  still  need  to  determine  V\  and  V"+1 .  We  want  V,  =  v*[h,  t„+i/a)  =  v*(0,  fn+1/4) 
and  so  the  appropriate  value  comes  from  the  second  equation  of 

3*(0,  <„+!/<)  g(tn+ 1/4)  +  \k[AfA~l  -  /)/(<„), 

i.e.,  # 

v  1  —  ffa(*n+j/4)  +  ifa-fcr,«^(2fagi(fn)  +  fie2ffa(*«))» 
where  g  =  (gi,ga)T-  Similarly, 

V'?+1  =  ga(<„+3/4)-  ?(2T^^7)(2£ag,1(tn+))  +  f,£2g/a(<n+,)). 

Computations  confirm  that  those  boundary  conditions  give  an  0{ik2)  globally  ac¬ 
curate  split  scheme.  Actually,  for  this  particular  example  with  k  —  2 h,  even  greater 
accuracy  can  be  achieved.  Computing  li ,(k)  from  (1.13),  the  truncation  error  for  Lax- 
Wcndroff,  shows  that  the  0(£fc3}  terms  exactly  cancel  the  0(tfc3)  terms  in  Z?aPiit(fc),  an<^ 
that  the  total  truncation  error  ETaM(k)u  is  actually  0(t2k3),  giving  0(eika)  global  ac¬ 
curacy.  Hy  retaining  more  terms  in  the  above  boundary  expansions  we  can  match  the 


Fig.  A  A.  Errors  in  the  computed  solution  U  of  Example  4.1  using  theO(ek a) 
boundary  conditions.  The  interior  error  is  0(cak2).  The  errors  arc  shown 
on  a  logarithmic  scale  For  various  values  of  k  and  c.  In  all  cases  t  =  0.2. 


FIG.  4.5.  Errors  in  the  computed  solution  U  or  Example  4.1  using  the  0{tk*) 
boundary  conditions.  The  interior  error  is  0(c*fca).  The  errors  arc  shown 
on  a  logarithmic  scale  hr  various  values  of  k  and  e.  In  all  cases  t  —  0.2. 


error  in  the  interior  solution.  Taking  one  more  term,  as  in  (4.27),  gives  0(ck3)  boundary 
data.  Figures  4.4  and  4.5  show  the  some  sample  results  using  0(eka)  and  0(ek3)  accurate 
boundary  data  respectively.  Errors  in  the  first  component  U  are  shown  at  time  t  =  0.2. 
Errors  resulting  from  the  boundary  conditions  have  propagated  in  to  approximately  x  — 
0.4. 

The  oscillations  in  the  error  near  the  boundary  are  due  to  the  fact  that  some  of 
the  boundary  conditions  used  (e.g.,  for  U0  )  have  sero  error  while  others  (e.g.,  for 

U0)  have  large  errors.  Since  the  split  scheme  is  only  mildly  dissipative  due  to  the  0(e) 
coefficients  in  the  Lax-Wendroff  step,  these  oscillations  introduced  at  the  boundary  die 
out  very  slowly  as  the  wave  propagates  into  the  interior.  This  is  in  no  way  an  indication 
of  instability.  Stability  for  this  example  follows  from  the  general  results  of  Section  4.0. 

4.3.  Variable  coefficient  systems — inflow  boundaries. 

Defining  the  proper  boundary  data  for  variable  coefficient  problems  is  not  significantly 
more  difficult  than  for  constant  coefficient  problems.  The  only  complication  comes  in 
switching  between  x-  and  ^-derivatives.  Consider  the  system  of  equations 

=  A(x,  t)  tt,  (4.29) 

and  for  simplicity  suppose  that  A/  is  constant,  while  A,  =  Aa(x,  t).  Proceeding  as  in 
(4.24), 

«*(0,  <„  4-  fc/2)  =  tt(0,  t„)  +  \kAfux(Q,  tn)  +  \kaA)u„( 0,  t„)  +  •  •  •. 

Now  we  must  be  more  careful  in  switching  back  to  t-derivatives.  We  have 

tt,(0,  tw)  =  A-‘(0,  tn)ut(0,  tn)  (4.30) 

and  by  differentiating  (4.29)  we  find  that 


t*tt  =  At  u*  +  Ati„ 

U(S  =  A*  u*  +  Au„ 

so  that 

«**  =  A-l(A_1(«tt  -  A,A-,ut)  -  A,A-1ut). 

Higher  order  derivatives  can  be  computed  similarly.  Continuing  as  in  (4.25),  we  obtain 

U*(0,tn)  =  ff(t„)+  ^A/A-l(0,tnV(fn)+  JkMjA-^O.fnHA-'^tn^fn) 

-  (A-'(0,  fn)At(0,  tn)  +  A*(0,  t„))A-'(0,  tn)</{tn)\  +  0{k»). 

This  can  be  truncated  in  the  usual  manner  to  obtain  an  appropriate  expression  for  U J. 
Example  4.3.  Consider  the  standard  quarter  plane  problem  for  the  scalar  equation 

ut  =  -(1  +  ent(x))us 
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with  Af  =  —I,  A,  =c  — ca(z)  and  c  <  1,  |a(z)|  <  1.  Then  (4.31)  gives 

.•(0. ,.) = *„) + l<rd*o)>,‘-) + 

x^U  +  ^WM  +  ot*1). 

We  thus  find  that  the  boundary  condition 

U*0  =  g(tn  +  k/2(l  +  «*(0))) 

is  0(tk 9)  accurate.  By  retaining  the  next  term  of  the  expansion  as  well  we  obtain  the 
0(c*3)  accurate  boundary  data 

vi = «((. + */2(i + »(o)» + i*a((,  ”!(!))>  V'1-)' 

The  other  necessary  boundary  data  can  be  generated  in  a  similar  manner. 


4.4.  Inflow-outflow  boundaries. 

Next  we  consider  a  constant  coefficient  problem  it*  =  Aux  Tor  x,  t  >  0  with  an 
inflow-outflow  boundary  at  x  =  0.  This  means  that  A  has  both  positive  and  negative 
eigenvalues.  For  simplicity  we  suppose  that  A  is  in  block  diagonal  form, 


0 

A11] 


(4.32) 


with  the  eigenvalues  of  A1  negative  and  those  of  A11  positive.  Partition  3  =  (u,  t>)T 
conformally  with  A.  Then  at  x  =  0  the  elements  of  u  are  inflow  variables  while  those  of 
v  are  outflow  variables.  The  boundary  conditions  are  assumed  to  be  of  the  form 


u(0,  t)  =  St/(0,  t)  +  g(t)  (4.33) 

where  5  is  a  constant  matrix  and  g  is  a  given  function.  We  now  split  A  as  A  —  Af  +  A, 
with  A/  and  A,  again  block  diagonal.  Moreover  we  suppose  that  the  eigenvalues  of  A/ 
arc  negative  and  those  of  A™  positive. 

* 

We  consider  only  the  problem  of  computing  U0  and  will  suppose  that  exp(JkA/dx)  is 
known  exactly.  Then  V0  *8  determined  from  the  interior  and  we  need  only  specify  Uq. 
As  usual,  we  introduce  3  (z,  t)  which  solves  the  subproblem  3t*  =  A/3*  and  find  as  in 
(4.24)  and  (4.25),  that 

3*(0,  tn  +  k/2)  =  3(0,  tn)  +  ikA/A-%1 0,  t«)  +  i*aAjA~a3«(0,  *,)  +  •  •  •. 

For  simplicity,  suppose  that  AaA-a  =  /  +  0(c).  Then  for  0((k *)  accurate  boundary 
conditions  we  can  take 

O'o  =  3(0,  tn+1/3)  +  ^(A/A-‘  -  /)3«(0,  tn).  (4.34) 
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Introducing  the  matrix 


B  =  A, A-  - 1  =  -  '  An{A,1r,  _  , 

we  can  rewrite  (4.34)  as 

Uo  =  «(0,  *n+l/a)  +  *n) 

V0  —  t>(0,  t«+i/i)  +  j^aawt(0,  <«). 

By  diiTerentiating  the  boundary  conditions  (4.33)  we  obtain 

“t(0,  *«)  =  5ut(0,  t„)  +  fl'(t„). 


Bn  0 
0  B„» 


(4.35a) 

(4.35b) 


Using  this  and  (4.33),  (4.35a)  becomes 

U0  —  (5v(0,<n+i/j{)  +  0(<n+i/a]  +  i^Bn[Svt(0,  tn)  +  ff'ft,,)].  (4.38) 


Recall  that  Vq  is  already  known.  We  can  thus  solve  (4.35b)  For  v(0,  tn+i/a)-  Using 
this  in  (4.36)  yields 


V o  =  S[V*0  -  \kBnvt  1  +  <7(t„+i/a)  +  **Bn|St>t(0,{n)  +  ^(t,,)] 
=  SV  „  +  ff(<n+i/a)  +  ik[Rutf(tn)  +  {BnS  —  SBn)vt(0,  tn)|. 

The  vt  term  must  in  general  be  approximated  by  a  finite  difference,  c.g., 

U0  =  SVQ  +  g(tn+i/a)  +  $kBug'(tn) 

+  (//,, 5 -55a2)(K5-Vr‘)  1- 


(4.37) 


(4.38) 


Alternatively  we  can  replace  vt  by  Alrvz  and  approximate  this  by  a  finite  difference  of  V 
at  time  tn.  This  approach  is  particularly  useful  when  more  terms  of  the  scries  are  kept 
and  higher  order  derivatives  must  be  approximated. 

The  use  of  such  boundary  conditions  is  illustrated  in  the  next  section,  where  the 
one-dimensional  shallow  water  equations  arc  considered. 

Boundary  data  at  points  near  the  boundary  can  be  found  in  a  similar  manner.  For 
example,  if  data  i/"+l  is  needed  for  some  0  <  j  <  p  we  can  expand  u  in  x-derivatives, 
switch  to  {-derivatives  along  the  boundary,  convert  these  to  {-derivatives  of  v  using  — 
Svt  +  gf ,  and  finally  switch  back  to  x-derivativea  of  v,  obtaining 


«(i»tn+ 1)  =  «(0,{B+i)  +  jh(Ar)  '(SA'^O.tn+i)  +  ,  . -x 

+  |ifc)V)-J[5(/l'f)Jr„(0(lB+1)  +  /(W,)]  +  "-.  1  ' 

These  boundary  conditions  are  suggested  by  Goldberg  &  Tadmor[21][22)  for  general 
inflow-outflow  problems. 


A 


The  shallow  water  equations  In  order  to  illustrate  the  derivation  of  intermediate 
boundary  conditions  for  a  more  realistic  example,  we  will  consider  the  one-dimensional 
shallow  water  equations  on  a  strip, 


.  •*•“•** 

with  initial  conditions 

2(*,0)  =  /(*),  0<x<l, 

and,  for  example,  the  boundary  conditions 

«(!,*)  =  0(l,t)-0o. 


(4.40) 


(4.41) 


(4.42a) 

(4.42b) 


Here  0o  is  the  mean  value  of  0  as  in  Section  2.0  and  the  boundary  condition  (4.42b) 
represents  nonrcflcction  at  the  boundary  x  =  1.  At  the  boundary  z  =  0  we  have  chosen 
to  prescribe  0.  Other  boundary  conditions  can  be  handled  similarly. 

As  in  Section  2.9,  the  equations  (4.40)  can  be  written  in  the  characteristic  form 
(2.46).  As  usual  we  suppose  that  |u|  <£  |0|.  Then  the  Ricmann  invariant  u  +  0  always 
flows  to  the  right  with  velocity  0/2  +  u  while  the  Riemann  invariant  u  —  0  always  flows 
to  the  left  with  velocity  0/2  —  u. 

The  problem  of  specifying  intermediate  boundary  conditions  is  simplified  if  we 
change  variables  and  compute  directly  in  terms  or  the  characteristic  variables,  which 
we  denote  by  p  and  a: 

p(x,  t)  =  u(x,  t)  +  0(x,  t), 

*(M)  =  «(*,*) -0(x,t). 

We  can  always  transform  back  to  find  u  =  [p  +  e)/2  and  0  =  (p  —  a)/2.  Rewriting  the 
differential  equation  (4.40)  in  terms  of  p  and  o  gives 


—  _  l[3P  +  1 

,  4l  o 


°  P 
p  +  3o  a  t 


(4.42) 


which  we  split  as  in  (2.48)  by  taking 


a.  —  if— ®  a _ i  3p  +  o  —  20o  0 

1  *[  0  0oJ*  J[  0  p  +  3o  +  20o / 


The  boundary  conditions  (4.42)  become 


p(0,  t)  —  o(0,  t)  +  2p(f) 

*(M)  =  -00- 


(4.43a) 

(4.43b) 


At  the  left  boundary  p  is  the  inflow  variable  and  the  boundary  condition  (4.43a)  is  of 
the  general  form  (4.33).  At  the  right  boundary  o  is  the  inflow  variable.  The  boundary 
condition  (4.43b)  indicates  that  a  is  constant,  and  hence  the  outgoing  wave  is  not 
reflected. 


i 
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For  k  =  4h/<j>0,  the  split  scheme  on  0  <  *  <  1  with  h  =  t/N  is  simply 

i?m  —  m  —  1)  2, . . . ,  N 

Sm~Sm+ 1»  m  ~  ~1»0, AT  —  1 

*  =LVF(A.,jfc)r|  ,  *n  =  0, 1, . . . ,  AT  —  1 

J?*+1  =  «m-l>  rn  =  1,2,...  ,N 

5»+‘=C+u  m  =  0, 1, . . . ,  TV  —  1. 

At  the  left  boundary  it  appears  that  we  need  to  specify  /?JJ,  iZlj,  /?J+1  and  Sg*.  In  fact 
So*  is  not  used  in  computing  SB+1  and  so  we  only  need  to  specify  the  R  values.  Note 
that  by  specifying  J?lt  we  avoid  having  to  specify  any  boundary  values  for  R" . 

The  given  boundary  conditions  (4.43a)  provide  Ro+i, 


RZ+i  =SS+1+2ff(fn+1). 


(4.43) 


Wc  next  apply  the  procedure  of  section  4.4  to  compute  Rq.  The  expression  (4.34)  provides 
0((k2)  accurate  boundary  data  for  the  quasilinear  problem  provided  A~l  is  evaluated  at 
(p(0,  t„),cr(0,  <„)).  The  matrix  B  —  A/A~l  —  I  is  given  by 


B  = 

and  the  expression  (4.37)  becomes 


0+0 

0 


-  i 


o 

-3*o  _  , 
p+io  1 


=  0(e) 


R°  ~  So  +  2®(i-+»/3)  +  **[((1^5??  3^)a|(0,  <n)  +  2(3p  +  o 

=  s'o + 2»(*» + + ff)) +  sfc((3p8+°il(p+Lj)'7‘(0,  tn) 


where  p  and  a  are  evaluated  at  (0,  t„),  This  can  be  approximated  by 

*o  =  So  +  2 g(tn  +  ak)  +  ~  $5" ‘) 

where 


(4.44) 


a  = 


3«3  +  S3' 

In  order  to  find  /?!,  wc  approximate  p*(—h,t„  +  kj 2).  This  is  equal  to  p*(0,t„  +  k) 
and  proceeding  as  in  Section  4.4  we  find  the  approximation 

R- 1  —  s\  +  2p(f«  +  2 ak)  + ^  +  3^n”^(5o  “  So  l)  (4.45) 

with  a  ns  above. 
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At  the  right  boundary  we  still  need  to  specify  Sq,  Sq*  and  So+I.  Since  the  boundary 
condition  (4.43b)  is  time-independent,  applying  the  general  procedure  at  this  boundary 
yields  simply 

So  =  So  =  S5+1  =  -*>•  (4.46) 

Figure  4.6  shows  the  result  of  some  computations  using  the  boundary  conditions 
(4.43),  (4.44),  (4.45)  and  (4.46).  The  boundary  conditions  at  the  right  boundary  have  not 
affected  the  interior  solution.  Errors  do  arise  at  the  left  boundary,  but  these  are  seen  to 
be  the  same  order  of  accuracy  as  the  interior  solution. 

As  in  Example  4.1,  the  oscillations  near  the  boundary  are  due  to  the  different 
boundary  conditions  being  of  different  accuracy. 

4.6.  Stability  of  the  initial- boundary  value  problem. 

In  general  stability  theory  for  initial-boundary  value  problems  is  considerably  more 
complicated  than  for  pure  initial  value  problems.  Only  recently  has  a  general  theory 
been  developed.  The  fundamental  paper  on  this  subject  is  by  Gustafsson,  Kreiss  & 
Sundslrom(30j. 

We  will  first  consider  an  inflow  boundary  with  boundary  conditions  as  derived  in 
Section  4.2.  In  this  situation  stability  can  be  proved  directly  from  the  Cauchy  stability  of 
the  interior  scheme  without  resorting  to  the  theory  of  Gustafsson,  Kreiss  and  Sundstrom. 
This  is  because  the  boundary  conditions  we  are  considering  are  independent  of  the  interior 
solution.  Consider,  for  example,  the  expression  (4.25).  Our  approximation  U q  can  be 
bounded  a  priori  in  terms  of  an  appropriate  discrete  Sobolev  norm  of  the  given  boundary 
data  </(f).  The  same  is  true  of  the  other  required  boundary  data. 

Stability  of  the  time-split  method  can  then  be  proved  using  the  following  general 
theorem,  which  states  that  any  Cauchy  stable  scheme  is  also  stable  for  the  initial¬ 
boundary  value  problem  provided  that  the  specified  boundary  data  {!/!),}„= o  >s  inde¬ 
pendent  of  the  interior  solution. 

THEOREM  4.1.  Suppose  Q(k)  is  Cauchy  stable.  For  the  initial-boundary  value 
problem  define  Un+i  by 

r/n+i  _  /QWm  rn>  p, 

\CS+l  m  =  0,l,...,p. 

Then  the  approximation  is  stable  in  the  sense  that 

IK'12  <  KT\\U°\\l  +  ATr||C||*  for  nk  <  T,  k  <  k0,  (4.47) 

where  Kt  and  K r  arc  constants  depending  only  on  T. 

Here  the  following  norms  are  used: 

\\in\l  =  h  f;  \irj\ 

m=0 
T/k  p 

l«  =  *EE  lG’la- 

,•=. 0 
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Proof.  By  the  Cauchy  stability  of  Q  (Stability  Definition  3.1'),  there  exists  a  constant 
a  and  a  norm  ||  •  ||,  equivalent  to  the  norm,  such  that 

ll<?{*)ll*  <  1  +  a*  for  k  <  ko.  (4.48) 

Extend  the  given  initial  data  {£/£»}  “=0  m  by  setting 

1/^  =  0,  m  =  — 1,— 2 . 


Then  solving  the  quarter  plane  problem  is  equivalent  to  solving  the  Cauchy  problem  and 
then  redefining  at  each  step.  Specifically,  we  set 


v«+i 


I'm  =Q(k)U 


m  =  0,  ±1,  ±2, ... 


and  then  take 


(/n+.  =  /G-‘  m”0'1 . "> 

otherwise. 


The  resulting  constitute  the  solution  of  the  quarter-plane  problem. 

By  (4.48)  we  have 

||i/n+1||a  <  (l  +  aA:)||r/n||8. 

By  (4.49)  we  obtain  the  following  bound  for  t/n+l: 

\\un+l\\ 2  <  ||t/”  M!i2  +  ||Gn+1||2 

where 

l|G"+l||a  =  AE|G7+1|*. 

i= 0 

Combining  (4.50)  and  (4.51)  gives 

||l/B+,||2  <  (1  +  «*)l|tH|2  +  !|Gn+,||2 

so  that  by  induction  we  obtain 


(4.49) 

(4.50) 

(4.51) 


||l/"||2  <  (l  +  a*n|l/°||2  +  £(1  +  afcn|GB-«||* 

9=0 

<«ar(l|G0!l2+Ellc—'||2) 

ror  nk  <  T.  Since  \\Un\\3  <  ||(/"||2,  ||f/0I(2  =  ||(/°||2  and 

£  l|Cf"-*||*  =  />££  |G*|2  <  £||C||? 

9=0  q^-Oj—O 

for  nk  <  T,  wc  obtain  the  desired  bound  (4.47)  with  K ?  —  eaT  and  Kt  ~  eaT/\.  | 
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To  sec  how  this  theorem  applies  to  a  time-split  method,  consider  the  method 


^n+1  = 

t/m+>  =  Q2{k)Um”+1- 


(4.52) 


We  have  added  the  index  n  +  1  to  U*  for  reasons  which  will  soon  be  apparent.  Suppose 
that  the  boundary  data  arc  of  the  form 


f/*n+1  =  G*n+I, 

f/”+1  =  G?+l, 


J  —  0, 1,  •  •  • ,  p, 
j  =  0,  l,...,p. 


(4.53) 


For  convenience  we  have  assumed  that  the  same  number  of  boundary  conditions  are 
needed  for  both  U'n+l  and  Un+l,  but  this  is  not  essential.  The  quanities  C*n+l  and 
G’"+l  are  determined  as  in  Section  4.2  in  terms  of  the  given  boundary  function  gr(<)  and 
some  of  its  derivatives  (say  d  derivatives).  Suppose  that  the  corresponding  Sobolev  norm 
of  </(<)  is  uniformly  bounded  by  some  constant  q: 


E  <  T 


Then  we  have 


l|C*ll?  <  Kn 


(4.54a) 


IlCIlf  <  K 27 

for  some  constants  K\  and  K 2. 

In  order  to  apply  Theorem  4.1  we  rewrite  (4.52)  as 


1  °]Mn+I  =  [°  <?«(*)' 

~Q*(k)  I  u  „  0  0 


(4.54b) 


ir  *in 

u 

.  U/  L 


(4.55) 


to  obtain  a  Cauchy  stable  scheme  for  the  “super-vector”  (U',U)T.  Note  that  the  method 
is  formally  implicit  oven  if  the  original  method  was  explicit,  as  it  must  be  since  the 
boundary  conditions  specified  for  (/*”*'  affect  the  computation  of  l/n+1.  The  Cauchy 
stability  of  (4.55)  follows  from  the  Cauchy  stability  of  Qv(k)Q \{k),  which  gives 


together  with 


IKHI  < 


||c ‘”11  <  c(||e°| 


where  C\  ~  C||Q|(fc)||.  Using  Theorem  4.1  and  the  bounds  (4.54)  we  find  that  (4.55)  is 
stable  for  the  initial-boundary  value  problem  and  that,  in  particular, 


<  /vT||f/°||2  +  kr[Ki  +  Ki)r 
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Inflow-outflow  problems.  Stability  can  also  be  demonstrated  for  inflow-outflow 
problems  with  boundary  conditions  of  the  form  discussed  in  Section  4.4.  As  above  the 
time-split  nature  of  the  scheme  can  be  handled  by  introducing  super-vectors.  Hence  we 
will  only  discuss  the  stability  of  a  general  one-step  scheme  in  which  the  inflow  variables  U 
and  the  outflow  variables  V  are  coupled  only  through  the  boundary  conditions.  As  usual 
we  assume  Cauchy  stability.  Our  discussion  will  be  rather  brief  but  similar  arguments 
can  be  found  in  Goldberg  &  Tadmor[21][22|. 

The  scheme  for  V  is  independent  of  U  and  we  will  assume,  as  we  did  in  Section 
4.4,  that  the  time-split  method  yields  a  one-sided  scheme  for  V  so  that  no  boundary 
conditions  need  be  specifled.  Then  from  Cauchy  stability  we  clearly  have 

\\vn\\l  <  \\v°\\i 

since  the  introduction  of  the  boundary  docs  not  affect  the  computation  of  {V”}?L0. 
Moreover  such  a  scheme  for  V  is  also  stable  in  the  sense  of  Definition  3.3  of  Gustafsson, 
Kreiss  &  Sundstrom[30]  (we  refer  to  this  as  GKS-stability).  This  stability  condition  also 
requires  bounds  on  a  norm  of  V  along  the  boundary.  The  GKS-stability  follows  easily 
from  the  theory  of  (30]  for  a  one-sided  scheme. 

GKS-stability  of  the  outflow  problem  is  just  what  we  need  to  prove  stability  of  the 
inflow  problem.  Recall  from  Section  4.4  that  the  boundary  conditions  for  U  depend  only 
on  g(t)  and  on  values  of  V  along  the  boundary,  and  can  be  bounded  in  terms  of  |||(7|||d 
and  || VU* .  The  former  of  these  is  assumed  to  be  uniformly  bounded  while  the  latter  is 
bounded  by  the  GKS-stability  of  V.  Theorem  4.1  thus  applies  to  the  inflow  problem  and 
hence  the  entire  approximation  is  stable  on  the  initial-boundary  value  problem. 

These  stability  results  arc  supported  by  large-time  numerical  calculations  for  all  of 
the  examples  which  have  been  given  in  this  chapter,  including  the  boundary  conditions 
of  Section  4.5  for  the  shallow  water  equations. 
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5.  Other  applications  of  the  theory 


5.1.  Introduction. 

In  this  chapter  some  of  the  theory  developed  in  previous  chapters  is  applied  to  a 
few  different  problems.  In  Section  5.2  hyperbolic  problems  in  two  space  dimensions  are 
considered.  Again  we  split  between  fast  and  slow  subproblcms  although  now  spatial 
splittings  may  also  be  used.  Intermediate  boundary  conditions  arc  derived  for  a  scalar 
example. 

We  then  turn  to  the  use  of  time-split  methods  for  problems  which  are  not  hyperbolic, 
since  many  of  the  techniques  that  have  been  introduced  are  applicable  to  other  problems 
as  well. 

In  Section  5.3  the  convection-diffusion  equation  ut  =  —  cux  +  cuxx  is  studied  as  a 
model  for  general  equations  containing  both  hyperbolic  and  parabolic  terms.  An  analysis 
very  similar  to  that  of  Section  2.5  is  performeo  with  analogous  results.  For  the  Cauchy 
problem  the  time-split  method  is  more  accurate  provided  the  mesh  ratio  is  chosen  ap¬ 
propriately.  For  boundary  value  problems  the  correct  intermediate  boundary  conditions 
at  the  inflow  boundary  can  be  computed  using  the  general  procedure  of  Chapter  4.  At  the 
outflow  boundary  no  special  boundary  data  need  be  specified,  but  the  solution  generally 
has  a  boundary  layer  at  this  boundary  which  causes  special  difficulties.  The  interior 
solution  (away  from  the  boundary  layer)  can  still  be  calculated  more  efficiently  than  with 
the  unsplit  method,  hut  less  efficiently  than  in  the  Cauchy  problem  due  to  mesh  ratio 
restrictions  imposed  by  the  boundary  layer. 

In  Section  5.4  a  very  different  kind  of  time-split  method  is  considered.  The  Pcaccman- 
Rachford  ADI  method  for  flic  two-dimensional  heat  equation  ut  —  uxx  +  uyy  is  viewed 
as  a  time-split  method  with  the  splitting  (1.43).  Ry  means  oT  the  procedure  of  Chapter  4, 
intermediate  boundary  conditions  arc  derived  for  a  rectangular  region  which  agree  with 
the  classical  boundary  conditions  for  this  method. 

5.2.  Hyperbolic  problems  in  two  space  dimensions 

The  time-split  method  can  be  used  in  two  (or  more)  space  dimensions  in  much  the 
same  way  as  in  one  dimension.  Locally  one-dimensional  methods,  where  a  time  . split 
method  is  used  to  reduce  a  multidimensional  problem  to  a  sequence  of  one-dimensional 
problems,  have  already  been  discussed  in  Chapter  I.  The  techniques  which  have  been 
developed  in  the  intervening  chapters  are  applicable  to  such  splittings  and  can  he  used  to 
analyze  their  efficiency  and,  in  some  cases,  to  generate  boundary  data  for  the  intermediate 
solutions.  This  will  be  done  for  the  I’eaceman-Raehrord  ADI  method  in  Section  5.4. 

In  this  section,  however,  wc  continue  to  concentrate  on  hyperbolic  problems  which 
can  be  split  into  “last"  and  “slow”  subproblcms.  Much  of  these  subproblcms  will,  in 
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general,  still  be  a  two-dimensional  problem.  In  some  eases  it  will  prove  useful  to  also  use 
spatial  splittings  in  order  to  solve  one  or  the  other  of  these  subproblems. 

A  general  hyperbolic  problem  in  two  space  dimensions  has  the  form 

=  Aiix  +  Buy  (5.1) 

where  the  matrices  A  and  B  have  the  property  that  £A  -1-  qB  is  diagonalizable  with 
real  eigenvalues  for  all  real  values  of  £  and  q.  In  the  notation  of  Chapter  1,  we  have 
A(u)  =  Aux  +  Buy  and  we  consider  splittings  of  the  form 

Ai(y)  =  AfUx  +  Bfuy 

yj2(«)  =  A,  u,  +  B,Uy. 

For  the  constant  coefficient  case  the  splitting  error  is  easily  computed  by  expanding  the 
exponential  solution  operators.  Define  the  differential  operators  C /  and  C ,  by 

Cf  =  Afdx  +  Bfdy 
.  C.  =  A,dx  +  Bxdy 


The  splitting  error  is  found  to  be 

=  -  s*3(K/C.  -  hc,c.c,  +  \c,c) 

~  \C\Cf  +  C'CfC.  -  $C/CJ)  +  0(k*).  1  '  ] 

which  is  analogous  to  the  one-dimensional  result  (2.22).  In  particular  the  splitting  error 
is  zero  if  all  of  the  matrices  Af,  Aa,  Bf  and  B,  commute. 

As  in  the  one-dimensional  case,  a  splitting  of  the  form  (5.2)  will  be  useful  if  A/  and 
Bj  are  sparse  relative  to  A  and  B  and  if  A„  and  B,  have  relatively  small  eigenvalues. 
Suppose  that  ||A/||  as  ||/f/||  aa  a  while  ||/ls||  aa  ||5,||  aa  to  with  e  <  1  and  that  the 
spectral  radius  of  each  matrix  is  comparable  to  its  norm. 

Lot  LW(A,  B,  k)  denote  the  two-dimensional  Lax-Wendroff  operator,  which  is  analo¬ 
gous  to  ( 1 . 1 1 )  and  can  be  found,  for  example,  in  Mitchell  and  Criffiths[41j.  The  stability 
limit  for  LW(A,  B,  k)  is  given  by 

r  ma \{p(A),p(B))  < 
n  \/S 

The  split  scheme  corresponding  to  (2.18a,c)  is  given  by 

Q,(k)  =  LW(A',B.,k) 


and 

Qf(k/2)  =  (LW(Af,Uf,k/m)r'2. 

An  efficiency  analysis  very  similar  to  that  performed  on  the  one-dimensional  problem  in 
Section  2.5  shows  that  the  optimal  mesh  ratio  on  the  fast  scale  is 


k  l 
mh  a 


(5-<) 
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This,  however,  violates  the  stability  limit  for  LW(Af,Df,k/m),  which  is 

jl<_l 

mh  -  V&a 

We  must  use  this  smaller  value  of  k/(mh)  instead.  Alternatively  we  can  replace  the 
two-dimensional  operator  LW(Af,Bf,k/m)  by  the  split  scheme 

LWx(Af,  k/2m)LWv(Bf,  k/m)LWx(Af,  fc/2m) 

where  LWX  and  LWy  represent  one-dimensional  Lax-WendrofF  in  the  x-  and  {/-directions 
respectively.  This  does  not  increase  the  truncation  error  significantly  but  increases  the 
stability  limit  to 

so  that  the  optimal  mesh  ratio  (5.4)  can  be  used.  (Recall  that  this  increase  in  the  stability 
limit  was  Strang’s  original  goal  in  introducing  the  Strang  splitting[49].) 

On  the  slow  scale  the  optimal  value  of  X  =  k/h  depends  on  the  size  of  the  splitting 
error.  If  J£apiu(k)  is  negligible  then 

X«  — . 

ea 

Again  this  violates  the  stability  limit  X  <  l/(\/8fa)  and  we  may  wish  to  introduce  a 
spatial  splitting  in  Q,(k). 

In  the  more  usual  situation,  however,  when  the  splitting  error  is  0(ea3fc3),  the 
optimal  mesh  ratio  is 

,  1 

X  - . 

y/ea 

This  is  well  within  the  stability  limit  and  there  is  no  need  to  introduce  spatial  splittings. 

Perturbed  problems.  The  splitting' (5.2)  is  also  useful  when  the  exact  solution 
operator  corresponding  to  >U(«)  is  known.  This  is  perhaps  not  so  common  in  two 
dimensions  as  in  one.  In  one  space  dimension  we  considered  several  examples  in  which 
the  coefficients  had  large  mean  values  and  small  variations.  We  could  then  pull  out 
a  constant  coefficient  subproblcm  ut  =  d/u*  which  could  be  solved  by  the  method 
of  characteristics.  Unfortunately,  in  two  dimensions  the  method  of  characteristics  is 
applicable  only  if  A/  and  Bf  arc  simultaneously  diagonalizable. 

Here,  however,  we  suppose  that  the  solution  operator  for  the  fast  part  is  known 
exactly  and  consider  the  time-split  method  (2.18a, b)  with 

Q.(fc)  =  LW(d„ /?.,*) 


QM  2)  =  cxp^k(Af<Jx  +  BfOy)). 

An  efficient  analysis  similar  to  that  of  Section  2.5  shows  that  the  optimal  mesh  ratio  is 
given  by 

X  =  —  if  the  splitting  error  is  negligible,  or 


X  =  — 
€a 

x  =  I 

a 


if  the  splitting  error  is  0(cfc3o3). 


In  the  former  case  we  can  only  acheive  Ihc  indicated  mesh  ratio  by  using  a  spatial  splitting 
for  Q«(&)  but  in  the  more  usual  situation,  when  splitting  errors  are  present,  this  is  not 
necessary. 

Boundary  conditions  for  a  perturbed  scalar  problem.  We  now  consider  a 
perturbed  scalar  problem  which  can  be  split  in  this  manner  and  show  how  to  derive 
appropriate  boundary  conditions  for  the  intermediate  solutions  in  two  space  dimensions. 

Consider  the  problem 

=  (1  +  *{x,  y,  t))us  +  (1  +  P{x,  y,  t))uv 
on  the  unit  square  (0,  lj  X  (0,  l)  with  boundary  conditions 

=  g\{y,t), 
u(x,l,t)  =  gu{x,t), 


(5.5) 

(5.6) 


and  suppose  that  |a(x,  y,  t)|  <  e,  \P{x,y,t)\  <  e  for  all  x,  y,  and  t  with  e  1.  It  is 
natural  to  split  this  by  taking 

Af  =  1,  A,  =  a[x,y,t), 

Bf  =  l,  BB=p(x,y,t). 


The  subproblem  ut  =  u*  +  tt*  can  be  solved  exactly: 

u*(z,y,  t  +  A:)  =  n'(x  +  k,y  +  k,t). 

Taking  k  =  2h  and  using  Lax-Wcndroff  on  the  slow  problem,  the  split  method  becomes 

Um,j  —  ^m+t.y+i  m,j  =  0, 1, . . .,  N  —  1, 

UZ,j  =  LW(a,  p,  k)U’mij  m,  j  =  1, 2, . . .,  N  -  1, 

UVJ  =  C+i.y+t  m,  j  =  0, 1, . . .,  N  -  1. 

The  values  U^J  and  arc  g*ven  by  the  boundary  conditions  (5.6)  while  the  values 

(7q*5  and  (/“0  are  not  required  in  computing  Un+l  and  therefore  do  not  need  to  be 
specified.  We  do  need  to  specify  the  following  intermediate  boundary  data: 

U  N,j>U  N,j  j  =  0,l,...,N, 

U  m,N>U  m,N  m  =  0,  1, . . .,  Af, 

First  consider  We  begin  as  usual  by  introducing  the  function  u*(z,y,  t) 

satisfying 

«t*  +  (5-7) 

with  initial  conditions  at  time  tn 

u‘(x,y,tn)  =  u[x,y,tn). 
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Then  we  want  U*Nj-  5=3  «*(1,  Vj,  +  */2).  Expanding  in  a  Taylor  aeries  and  UBing  (5.7) 
together  with  the  fact  that  a*  =  u*  and  u*  =  uy  at  time  tn,  we  obtain 

«*(l,y,{«  +  */2)  =  u*(l,  y,tn)+  kkut(l,y,tn)  +  |*aut*{(l,y,{„)  +  --- 

=  «(1  ,y,t„)+  5*(«x(l,  V.  tn)  +  uy(l,  y,  *„))  (5-8) 

+  i*2(«x«(l » y,  *«)  +  2uItf(l,  y,  tn)  +  uyv(l,  y,  <«))  +  •  •  •. 

We  now  use  the  original  equation  (5.5)  to  replace  z-dcrivatives  by  y-  and  {-derivatives, 
ao  that  the  given  boundary  data  gi(y,t)  and  its  derivatives  can  be  used  to  specify  U us. 
We  find  that 

,  fc  fc8  ( 

«  (t,y,<«  +  k/2)  =  u+  2(f+ Q)(u«  +  (®  -  0K>  +  +  ^u‘» 

+  (1  +  p)3uvy  +  +  Q*)u‘  +  "  °*) 

+  (1  +  P)(Py  —  «*)  —  Pt  ~  (i  +  +  •  •  * 

where  the  functions  u,  a,  and  p  (and  their  derivatives)  on  the  right  hand  side  are  all 
evaluated  at  (l,y,  t„).  If  a,  0  and  their  first  derivatives  arc  all  0(e),  then  this  can  be 
simplified  in  the  usual  manner: 

«*(1  ,V,tn  +  k/2)  =  u(l,y,<„  +  k/ 2(1  +  a))  +  >»><«)  +  0(cfca). 

We  can  maintain  the  accuracy  of  the  interior  scheme  by  using  the  boundary  values 

u*N,j  =  <?i(iMn  +  \k/(  1  +  a(l,jh,  {„)) 

f  k(cc(l,jh,tn)  —  P(\,jh,tn))\  t  s 

*  { - 2fTM7 

The  boundary  values  U ^  N  along  the  boundary  y  =  1  arc  found  in  exactly  the  same 
manner.  We  obtain 

Urn.N  =  32  (j'M«  +  ik/( !  +  0[jh,  1  ,{„)) 

,  ( k{fl{jh,  l,tn)-a(jh,  1,{„))V  ,,L  .  * 

To  compute  boundary  values  for  the  second  intermediate  solution  f/**,  we  proceed 
as  we  did  in  the  one  dimensional  case  by  defining  u  ( x,y,t )  as  the  solution  of 

**  **  *• 

Ut  =  UX  +  Uy 

with  u**(z,y,{n+i)  =  u(x,y,tn+i)  and  then  solving  backwards  in  time  to  find  U £ j  s=» 
u**(l, jh, tn+i  -  k/2).  The  expression  we  obtain  for  «**(!, jfi,tn|i  —  k/2)  is  exactly  the 


same  as  the  right  hand  side  of  (5.8)  but  with  A  replaced  by  —k  an{f  all  functions  evaluated 
at  (l,y,tn+i)-  We  can  thus  take 


Unj  =  0i(M«+i  -  $*/(!  +  ar(l,iM»+t)) 


fk[a(l,jh,tn+i)-0(i,jh,tn+i))\  „ 

V  2(1  +  at(l,jh,tn+i)) 


tjhy  l»+l) 


with  a  similar  expression  for 

Irregular  regions.  Attempting  to  compute  on  irregular  (nonrcctangular)  regions 
generally  complicates  the  problem  of  specifying  boundary  conditions  for  any  numeri¬ 
cal  method.  Gridpoints  frequently  do  not  lie  exactly  on  the  boundary  and  so  special 
procedures  must  be  used  for  points  near  the  boundary  even  when  the  correct  data  are 
known  along  the  boundary  itself.  Here  we  will  only  consider  the  problem  of  transforming 
boundary  conditions  for  the  given  problem  into  boundary  data  for  the  intermediate  solu¬ 
tions.  The  problem  of  then  using  these  data,  defined  along  the  boundary,  to  specify  the 
necessary  solution  values  at  nearby  points  can  then  be  handled  by  standard  techniques. 

Consider  the  problem  (5.5)  in  a  region  with  boundary  parametrized  by  (x(s),  y(s)), 
0  <  s  <  1.  The  region  is  assumed  to  lie  to  the  left  of  this  curve.  The  boundary 
conditions  arc 

u(x(s),y(a),t)  =  g(s,f) 

at  inflow  points.  For  convenience  we  will  assume  that  a  and  P  are  independent  of  t  and 
will  write  a(s)  for  a(x(s),  y(s))  and  similarly  for  /?.  Then  (x(a),  y(s))  is  an  inflow  point  if 


x'(s)(l  +/?(«))  <  y'(s)(l  +  «(«)). 


This  is  illustrated  in  Figure  5.1. 

For  the  rectangular  region,  we  replaced  x-derivatives  by  y-  and  t-derivatives  at 
the  right  boundary  while  at  the  top  boundary  we  replaced  y-dcrivativcs  by  x-  and  t- 
derivatives.  These  are  both  special  cases  of  the  general  situation.  At  the  inflow  boundary 
of  a  nonrcctangular  region  we  must  replace  both  the  x-  and  the  y-dcrivativcs  by  tangential 
and  ^-derivatives  in  order  to  obtain  expressions  in  terms  of  the  given  boundary  data. 

In  determining  u  at  inflow  points  we  first  obtain  an  expression  analogous  to  (5.8), 

u*(x(s),  y(»),  tn  +  A/2)  =  u(x(s),  y(s),  <n) 

+  £*(«x(z(s),  y(»),  tn)  +  «v(z(»),  y(«),  tn))  +  •  ”• 

Now  we  solve  for  ux  and  uv  in  terms  of  the  given  boundary  conditions  from  the  equations 

fftMn)  =  (*(*),  !/(»),<«) 

=  (l  +«(»))ttx(x(»),y(«),  <„)  +  (!  +/?(8))uy(x(s),y(»),<n) 

0.(Mn)  =  x\a)uz(x(s),y(s),tn)  +  y'(s)«„(x(s),y(  »),<„). 

Solving  this  system  gives  (dropping  arguments  for  clarity) 

_ _ y'(h  -  (1  +  P)g a 

U *  (1  +  P)x'  -  (1  +  a)y/ 

_  z'gt  -  (I  +  ot)g, 

Uy  (1  +  p)x'  -  (1  +  a)y' 
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FlG.  5.1.  Irregular  region  with  boundary  parametrised  by  (x(8),y(s)),  mov¬ 
ing  counterclockwise  as  8  increases.  The  characteristics  in  the  x-y  plane  are 
also  shown.  The  exact  solution  propagates  along  these  lines,  which  have 
slope  (1  4-  fi(x,y))/(l  +  a(x,y))  =  i  +  0(e)  at  each  point  ( x,y ).  The  inBow 
portion  of  the  boundary,  where  the  boundary  conditions  are  specified,  is 
shown  as  a  bold  line. 


so  that 


«*  4-  Uy 


(x/  -  y')gt  +  {fl~  a)ga 

(1  +  f))x'  -  (l  +  a)y' 


Note  that  the  denominator  is  nonzero  at  inflow  points.  Similarly,  we  can  obtain  expres¬ 
sions  For  higher  derivatives.  When  a,  and  their  derivatives  arc  ()(t)  we  have 


«*(*(»).  »(«)»*"  +  fc/2)  =  »(».  <*(«)) 


fc(/3(s)-a(s)> 


(I  +  0(s))x'(s)  -  (1  +  a(s))y'(s) 


where 


'*(«)  =  <«  +  ^( 


k{x'(8)  -  i/'(s)) 


(1  +  P(s))x'(s)  -  (1  +  n(»))y'(s) 


> 


This  Formula  reduces  to  those  derived  beForc  on  the  unit  square,  in  which  case  cither 
z^s)  =  0  or  y'(s)  =  0. 


5.3.  Convection-diffuaion  equations. 

As  mentioned  in  Section  1.6,  lime-split  methods  are  Frequently  used  to  solve  equa¬ 
tions  oF  mixed  type  by  splitting  between  the  hyperbolic  and  parabolic  parts  and  using 
dilFercnt  methods  on  the  two  pieces.  This  is  done  for  example  with  the  Navicr-Stokes 
equations  For  viscous  fluid  flow  or  convection-diffusion  equations  for  miscible  flow. 

In  this  section  we  consider  a  simple  model  equation  for  such  problems,  the  constant- 
coefficient  scalar  convection-diffusion  equation 


ut  =  —(Mix  +  tu„ 


(5.9) 
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where  c  and  r  arc  nonnegativc  constants.  Consider  the  splitting  A i(«)  =  —  cux,  yfa(ti)  = 
tuxx.  For  this  scalar  constant  coefficient  problem  there  is  no  splitting  error  so  we  do  not 
need  to  use  the  Strang  splitting.  If  k/h  —  —  p/c  for  some  positive  integer  p  then  the 
subproblem  uj  =  —  cu*  can  be  solved  exactly.  Using  Crank-Nicolson  for  the  remaining 
subproblcm  gives  the  split  method 

ul  = 

Ul+l=CN(e,k)U*m 
where  CN  is  the  Crank-Nicolson  operator,  which  has 


(5.10) 

the  form 


CN(A,  *)  =  (/-  ^kAD+D~)~l(I  +  \kAD+D~) 

for  a  general  constant  coefficient  system  ut  =  Auxx. 

If  we  eliminate  the  intermediate  solution  U* ,  we  can  rewrite  the  split  method  (5.10) 
as  a  one-step  method: 


(/  -  'skcD+D-)U^1  =  (/  +  \keD+D-)U*m_p.  (5.11) 

Figure  5.2  shows  the  stencil  for  this  method  when  c  —  p  —  1. 


The  method  can  thus  be  viewed  as  a  “skewed  Crank-Nicolson”  method  similar  to 
the  skewed  Lax-WcndrolT  method  of  Example  1.2.  The  stencil  of  the  method  follows  the 
characteristic  of  the  hyperbolic  part  of  the  problem. 

Time-split  methods  similar  to  (5.10)  can  be  used  for  more  general  systems  of  the 

form 

ut  =  Aux  +  Buxx 

where  A  and  B  may  be  functions  of  x,  t,  and  u.  One  way  to  proceed  is  to  use  the 
splitting  A i («)  =  Aux  and  A-i(u)  =  Buxx.  In  general  neither  subproblcm  can  be  solved 
exactly,  but  it  may  be  advantageous  to  use  different  numerical  procedures  for  the  two 
subproblems.  This  is  the  approach  generally  taken  with  the  Navicr-Stokcs  equations(l). 

Another  alternative  is  to  use  the  splitting  >?i(m)  =  A/ux  and  A^u)  =  Aaux  +  Buxx 
where  ut  =  A/ux  can  be  solved  exactly  and  the  remaining  subproblcm  ut  =  ^(u) 
corresponds  to  small  perturbations,  i.c.,  p(Aa)  and  p(B)  arc  small  compared  to  p{Af). 

Here  we  consider  only  the  scalar  problem  (5.9),  which  is  sufficient  to  illustrate  some 
of  the  new  issues  that  arise  when  time-split  methods  are  applied  to  such  problems.  In 
particular,  when  r  is  small  there  may  be  a  boundary  layer  at  the  outflow  boundary,  which 
poses  special  problems  for  the  time-split  method. 
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Efficiency  analysis  for  the  Cauchy  problem.  Before  considering  boundary 
value  problems,  however,  we  first  perform  an  efficiency  analysis  for  the  Cauchy  problem 
similar  to  that  of  Section  2.5.  The  split  method  (5.10)  will  be  compared  to  the  unsplit 
method 


(l  +  lkcDo-lkeD+D-py-' 

=  (1  -  ±kcD0  +  ^keD+D-) U^. 


(5.12) 


This  method  is  second  order  accurate  with  a  truncation  error 


E(k)u  —  fc{fc8[—  ^c3  +  jc2edx  —  jet2#2  +  ^e3#3] 
+  >»8[-  +  iVed*I}«***  +  0(k4). 


(5.13) 


We  ast  me  that  c  «  1,  (  <<  1  and  that  u  is  smooth  so  that  all  deri  vatives  of  u  are 
order  unity.  (This  latter  assumption  will  not  be  valid  in  the  boundary  value  problems 
considered  later.)  Then  E(k)  can  be  approximated  as 


E{k)  «  -  &k[k2c3  +  2h2c\d3s. 

We  sec  that  the  error  is  dominated  by  terms  arising  from  the  hyperbolic  part  of  the 
equation.  The  global  error  at  time  t  =  1  is  roughly  bounded  by 


i\\E(k)u\\  «  &*2c3  +  2ft2c)||uxxx||.  (5 .14) 

As  in  Section  2.5  the  optimal  mesh  ratio  is  found  by  requiring  Jfc2c3  «=s  2 h2c,  for  otherwise 
we  could  increase  one  of  the  stepsizes  without  significantly  increasing  the  error.  This 
gives  the  optimal  value  of  the  mesh  ratio: 


X  rw  \/2/c. 


(5.15) 


For  comparison  purposes  we  wish  to  normalize  the  error  at  t  =  1  by  the  amount  of  work 
performed.  A  tridiagonal  system  of  equations  must  be  solved  at  each  step  but  this  only 
requires  work  proportional  to  l / h.-  The  work  required  to  compute  the  solution  at  t  =  1  is 
thus  proportional  to  1  fkh.  The  same  is  true  in  the  split  method,  with  a  similar  constant 
of  propo  tionality,  and  so  we  can  normalize  the  error  simply  by  dividing  by  kh.  Using 
(5.14)  an  1  (5.15),  we  find  that 


normalized  error  «  ^c2||uxxx||.  (5.16) 

Now  consider  the  split  method  (5.10).  Since  there  is  no  splitting  error  and  no  error 
is  committed  in  solving  the  first  subproblem,  the  overall  truncation  error  is  simply  the 
truncation  error  of  CN(c,k),  which  is  found  by  selling  c  =  0  in  (5.13): 


ETS“{k)  u  =  i V*(*Vt?3  +  h2fds]uxxx  +  (){k*). 
The  optimal  mesh  ratio  is  thus 

c  V 


(5.17) 


(5.18) 


86 


This  indicates  that  large  timesteps  arc  optima!,  X  =  0(  1/c). 

Using  (5.18)  in  (5.17)  gives  the  following  normalized  error  at  t  =  1: 

normalized  error  »  g£2>/[|3£u||  ll^xull>  (5.19) 

This  is  smaller.than  (5.16)  by  roughly  a  factor  of  (c/c)2. 

These  results  are  virtually  identical  to  those  of  Section  2.5  for  the  pure  hyperbolic 
problem  in  the  same  situation,  namely  for  a  perturbed  problem  with  no  splitting  error. 
Other  situations,  e.g.  splittings  with  error  or  the  use  of  several  steps  of  a  difference 
method  on  the  fast  problem,  can  be  investigated  in  the  same  manner  with  analogous 
results.  Numerical  experiments  have  confirmed  these  theoretical  predictions. 

Nonsmooth  solutions.  The  advantages  of  the  time-split  method  are  most  clearly 
seen  when  computing  nearly-discontinuous  solutions,  for  example  shocks  in  the  Navier- 
Stokes  equations  or  steep  concentration  gradients  in  miscible  flow  problems. 

Example  5.1.  Consider  the  model  problem 


ut  =  —ux  +  euxx 


(5.20) 


with  initial  conditions 


x  <  0.1, 
x  >  0.1. 


This  initial  discontinuity  smears  out  as  it  propagates  to  the  right  with  speed  1.  At  time 
<  =  0.7  the  true  solution  is  seen  as  the  dashed  line  in  Figure  5.3  (with  e  —  10~3). 

The  unsplit  method  (5.12)  performs  poorly  on  such  problems  because  of  the  convec¬ 
tive  term.  The  resulting  solution  is  oscillatory  as  seen  in  Figure  5.3,  which  shows  the 
solution  obtained  with  k  =  h  —  10-2. 

With  the  split  method  (5.10)  the  convection  is  handled  exactly  by  shifting.  Only  the 
diffusion  is  handled  numerically  and  discontinuous  initial  data  cause  no  problems.  By  the 
efficiency  analysis  it  is  optimal  to  take  X  =  0(  1/t).  We  choose  to  again  take  h  =  10-2 
and  take  k  =  0.7  which  corresponds  to  X  =  70.  Figure  5A  shows  the  resulting  solution 
obtained  witii  a  single  step  of  the  time-split  method  (5.10). 


Boundary  value  problems.  We  now  turn  to  the  most  interesting  case:  the 
boundary  value  problem  (5.9)  on  0  <  x  <  1.  For  definiteness  we  will  take  c  =  1.  When 
f  =  0  the  equation  is  the  familiar  hyperbolic  equation  ut  =  -u,  for  which  boundary 
data  need  only  be  specified  at  the  inflow  boundary  x  =  0.  The  exact  solution  is  a  wave 
moving  to  the  right,  unaltered,  with  speed  1.  When  c  is  small  the  solution  is  again  a  wave 
which  moves  to  the  right,  but  now  it  dissipates  slowly  as  it  moves  along.  For  t  very  small 
we  might  expect  the  solution  to  be  very  similar  to  that  obtained  with  f  =  0.  However, 
whenever  <  >  0  the  equation  (5.9)  is  parabolic  and  boundary  conditions  must  be  specified 
both  at  x  =  0  and  at  x  —  1.  Equation  (5.9)  is  a  singular  perturbation  equation  since  the 
limiting  equation  with  c  =  0  has  a  singular  nature  quite  different  from  that  with  r  >  0. 

For  small  c  the  solution  to  (5.9)  has  a  boundary  layer  near  x  =  1,  a  small  region  in 
which  the  solution  changes  rapidly  in  order  to  match  the  boundary  conditions  at  x  =  1. 
For  (5.9)  the  boundary  layer  has  width  ()(f).  For  x  <  1  —  t  the  solution  is  simply  a 
rightward  moving  wave,  slowly  dissipating,  and  looks  very  much  like  the  soulion  to  the 
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FIG.  5.3.  Solution  of  the  convection-diffusion  equation  or  Example  5.1  with  (  =  10~3 
at  time  t  =  0.7.  The  dashed  line  is  the  true  solution.  The  solid  line  is  the  solution 
computed  with  the  unsplit  method  (5.12)  with  k  —  h  =  10~a. 


Fig.  5.4.  Solution  of  the  convection-diffusion  equation  of  Example  5.1  with  c  =  10~s 
at  time  t  —  0.7.  The  dashed  line  is  the  true  solution.  The  solid  line  is  the  solution 
computed  with  the  split  method  (5.10)  with  h  =  10-8,  k  =  0.7. 


FlG.  5.5.  True  solution  (dashed  line)  and  computed  solutions  to  the  steady- 
state  convection-diffusion  solution  using  the  time-split  method  (5.10)  with 
h  =  10~2  and  several  different  values  of  k.  Again  e  =  5  X  10-2. 


pure  hyperbolic  problem.  The  solution  is  smooth  and  the  convective  term  -u,  dominates 
the  dissipative  term  (uxx.  In  the  region  1  —  e  <  x  <  1,  however,  the  solution  is  rapidly 
changing  and  uxx  jiix.  Here  both  terms  arc  equally  important  and  the  solution  in 
this  region  is  quite  unlike  that  seen  for  t  —  0. 

The  presence  of  the  boundary  layer  in  this  problem  causes  difficulties  in  the  applica¬ 
tion  of  any  numerical  procedure.  The  time-split  method  performs  quite  well,  particularly 
in  the  interior,  provided  that  moderate  values  of  the  timestep  k  arc  used.  Using  X  i=»  1 
gives  results  which  arc  everywhere  better  than  the  unsplit  method,  by  a  Tactor  of  e  in 
the  interior  (see  Example  5.2).  This  is  to  be  expected  based  on  the  previous  analysis  of 
the  Cauchy  problem.  However,  it  is  no  longer  possible  to  obtain  even  greater  efficiency 
by  using  larger  tirnesl.eps.  This  is  because  or  errors  arising  in  the  boundary  layer.  It  is 
illuminating  to  analyze  the  difficulties  which  rise  when  larger  timesteps  are  used. 

In  order  to  concentrate  on  the  boundary  layer  itself,  we  lirst  consider  the.  steady-state 
problem  (5.20)  with  time-independent  boundary  conditions 

u(0,  t)  -  1  -  e~l/‘ 

«(M)  =  0 

and  initial  conditions 

u(z,0)  =  l 

The  solution  to  this  problem  is  simply 

u(x,t)  —  1  (5.21) 

for  all  t,  as  shown  by  I, lie  dashed  line  in  Figure  5.5  for  i  =  5  X  10— 2. 
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The  numerical  solution  to  this  problem  will  also  reach  a  steady  state,  though  in 
general  the  numerical  steady  state  will  be  different  from  the  true  solution.  For  the  unsplit 
method  (5.12)  setting  U JJ,  =  U\ J+*  shows  that  the  steady-state  solution  satisfies 

(7C0O  +  =  0. 

This  solution  depends  on  h  but  is  independent  of  the  timestep  k  used  to  compute  it  (and 
hence  is  independent  of  X).  The  numerical  solution  is  quite  accurate  if  h  is  sufficiently 
small.  If  h  >  £  then  oscillations  begin  to  appear  in  the  boundary  layer.  This  will  be  seen 
in  Example  5.2. 

Now  consider  the  time-split  method  (5.10).  In  order  to  implement  this  method  we 
must  specify  some  additional  boundary  values  Uq,  U t ...t/p_,  at  the  boundary  z  =  0. 
This  can  be  done  using  the  general  procedure  of  Chapter  4,  as  will  be  seen  later  in 
this  section.  For  the  present  example  with  time-independent  boundary  conditions,  these 
simply  reduce  to 

17*  =  1  —  e~1/<e  for  j  =  0,  t,...,p—  1. 

At  the  boundary  x  =  t,  where  the  boundary  layer  is,  we  do  not  need  to  specify  any 
additional  boundary  data.  Figure  5.5  shows  the  numerical  steady-state  solution  for  the 
time-split  method  with  <  =  5  X  10-2,  h  =  10-2  and  several  different  values  of  k.  Note 
that  this  steady-state  is  no  longer  independent  of  k  and  is  far  from  the  true  solution  even 
for  moderate  values  of  X. 

In  order  to  understand  this  phenomenon  we  must  consider  the  individual  steps  of 
the  time-split  method  in  more  detail.  In  the  first  step  of  (5.10)  the  solution  shifts  to  the 
right  and  much  of  the  boundary  layer  is  lost.  If  k  >  6  the  boundary  layer  shifts  almost 
completely  out  of  the  interval.  We  then  have 

«  1  (5.22) 

for  all  m.  The  solution  {/*  is  nearly  independr  ./  k  for  k  >  c.  In  the  second  step  of 
(5.10)  we  arc  using  Crank-Nicolson  to  solve  the  heat  equation  with  initial  values  (5.22) 
and  boundary  values  (/q  +  1  ~  1,  L/Jy"1  =  0.  For  large  values  of  k  the  solution  Un+X 
tends  to  t/£,+  l  =  1  —  zm,  which  is  the  steady-slate  solution  of  ut  =  <uxx  with  boundary 
conditions 

u(0,  t )  =  1,  u(l,  t)  =  0. 

This  explains  the  “over-diffused”  nature  of  the  solutions  seen  in  Figure  5.5  for  larger 
values  of  k. 

Where  does  the  time-split  method  break  down?  Recall  that  there  is  no  splitting  error 
for  this  problem  so  that  if  both  subproblems  arc  solved  exactly  we  should  obtain  the 
exact  solution  to  the  original  problem,  for  any  value  of  k.  Let  us  now  do  this.  The  first 
subproblem  ut  =  —  ux  is  already  being  solved  exactly  (modulo  the  boundary  conditions 
at  x  =  0,  but  these  have  a  negligible  effect  on  the  results  seen  here).  We  now  wish  to 
also  use  the  exact  solution  operator  for  the  second  subproblem 

u**(z,t)  =  fu**(x,t),  0  <  x  <  1,  t„  <  t  <  tn+ 1 

with  initial  conditions 

«“(x,  In)  —  u[x,tn{  i). 
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FIG.  5.6.  The  correct  boundary  conditions  u**(l,  t„  +  r)  for  0  <  r  <  k  with 
£  =  5  X  I0“2  and  k  =  1. 


In  order  to  apply  the  exact  solution  operator  we  must  also  specify  boundary  conditions 
u**(t,tn  4-  r)  for  all  r  with  0  <  t  <  k.  This  can  be  done  in  the  standard  way.  We  will 
use  the  fact  that  wc  know  the  true  solution  u(x,  <n+J)  =  u*{x,tn+\)  and  expand  about 

+  t)  =  w**(l ,<„+i)  +  (r  -  £)«“(! ,<„+i)  +  Kt  ~  k)2u*t*{l,tn+i)  +  ••• 

=  »**(!,  tn+i)  +  (r  -  *)f»**(Mn+i)  +  ^(r  -  k)2c2u*xxxx(\,tn+i)  +  •••. 

This  can  be  evaluated  directly  using  the  steady-state  solution  (5.21).  We  find  that  the 
proper  boundary  conditions  are 

u**(Mn  +  r)  =  l-c(T-fc>/'. 

This  is  shown  Tor  k  =  1  in  Figure  5.6.  Note  that  me  boundary  value  remains  nearly 
constant  through  most  of  the  time  interval.  Only  at  the  end  docs  it  suddetdy  drop  to 
the  final  value  u  (l,<n+t)  =  0.  This  explains  why  the  resulting  solution  is  not  “over- 
dilTuscd”  when  the  true  solution  operator  is  used.  No  diffusion  occurs  until  near  the  end 
of  the  time  interval,  for  k  —  c  <  r  <  k,  and  the  length  of  this  interval  is  independent  or 
k,  as  it  must  be  since  the  resulting  true  steady-stale  solution  is  independent  of  k. 

This  shows  where  the  time-split  method  breaks  down  for  large  k.  When  Crank- 
Nicolson  is  applied  in  the  second  step  the  correct  boundary  values  are  used  at  the  ends  of 
the  time  interval,  but  no  account  has  been  taken  of  the  nonsmoolh  behavior  of  u  (l,tn4- 
t)  for  0  <  r  <  k.  Since  Crank- Nicolson  is  only  accurate  for  smooth  boundary  data,  we 
get  poor  results.  It  is  as  if  we  had  used  the  exact  solution  operator  with  smooth  boundary 
data  obtained  by  linearly  interpolating  between  u  (1,<„)  =  1  and  u  (l,fn+i)  —  0. 

It  appears  that  this  difficulty  with  the  time-split  method  can  be  avoided  oidy  by 
taking  k  sufficiently  small.  If  k  <  r  then  the  boundary  layer  is  not  entirely  shifted  out 
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of  the  interval  and  the  resulting  true  boundary  conditions  for  «**  are  much  smoother. 
Moreover,  reexamining  (5.18)  shows  that  we  also  want  k/h  smaller  than  was  optimal  for 
the  Cauchy  problem.  From  the  steady-state  solution  (5.21)  we  see  that 

l|3j«ll » 1  Ab¬ 

using  this  in  (5.18)  gives  an  optimal  mesh  ratio  for  the  split  method  of 

X  sa  1 


rather  than  0(i/e)  as  was  optimal  for  the  Cauchy  problem. 

To  summarise,  we  find  that  for  boundary  layer  calculations  we  should  take  Jfc  s=r  h  < 
e,  for  the  split  method,  just  as  we  would  for  the  unsplit  method.  The  resulting  solution 
is  more  accurate  than  that  obtained  with  the  unsplit  method,  by  a  factor  of  e  in  the 
interior. 

Intermediate  boundary  data  at  the  inflow  boundary.  Now  consider  the 
unsteady  boundary  value  problem  (5.20)  with  boundary  conditions 


u(0,  t)  —  g[t) 

u(M)  =  0. 


(5.23) 


We  must  specify  additional  boundary  values  U'Q,  U\, .  Consider  the  step  from 

to  tn+ 1  and  let  u*  satisfy 


ut  —  -u*,  t  >  tn, 

u*[x,tn)  —  «(M»),  o  <  x  <  1. 

Then  we  want 

U)  =  u(jh,  tn  +  k) 

=  u*(0,  tn  +  k-  jh). 

Let  t  =  k  —  jh.  Expanding  about  u*(0,  t„)  and  using  (5.24), 

U ]  =  u*(0,  t„)  +  rut*  (0,  tn)  +  $7-au*t(0,  <n)  +  0{k3) 
=  «*(0,  <»)  -  tu* (0,  tn)  +  ^t2u*x(0,  t„)  +  0[k3) 
—  «(0,  tn)  -  rut( 0,  t„ )  +  %t2uxx( 0,  t„)  +  0(fc3). 


(5.24) 


(5.25) 


In  order  to  use  the  given  boundary  data  (5.23)  we  wish  to  replace  z-dcrivativcs  of  u  by 
{-derivatives  using  (5.20).  After  some  manipulations  we  find  that 

««  =  -  2£u,„  +  £au,„« 


and  so 


U*  =  -«t  +  cuxx 

=  -Ut  +  +  2eau„,  -  e3u«„ 

=  -ut  +  iutt  -  2tauttt  +  0{f  *). 
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Continuing  to  replace  z-dcrivatives  by  t-derivaUvea,  we  obtain  a  power  series  in  e  which 
involves  only  time  derivatives.  Similarly  we  find  that 


2et»m  +  0(€*). 

Using  these  expressions  in  (5.25)  gives  the  desired  expansion: 


Uj  —  t*(0,  tn)  —  t(— u*  +  cut*  —  2eauttt  +  •*•)  +  —  2cttttt  +  •  •  •) 

=  g( 0,  *«  +  t)  -  reg"{tn)  +  (2rca  -  ra«)pw,(tn)  +  0(e3k  +  *V). 

Note  that  this  approach  works  only  when  e  <  I. 

Example  5.2.  Consider  (5.20)  with  e  =  10-s,  initial  conditions 


(5.20) 


u(z,0)  =  1  -  z, 


and  boundary  conditions 

u(0,  t)  =  g(t)  =  1  +  0.1  sin(2  xt) 

«(1,  t)  —  0. 

Figure  5.7  shows  the  solution  at  time  t  =  2  using  the  unsplit  method  with  k  —  h  —  10-a. 
Figure  5.8  shows  the  results  with  the  time-split  method  (5.10)  again  with  J b  =  h  =  10-a 
(p  =  1)  using  (5.26)  for  U*Q.  Note  that  the  oscillations  in  the  boundary  layer  have 
disappeared.  Moreover  the  interior  solution  is  more  accurate  by  a  factor  of  e,  as  can  be 
seen  in  the  accompanying  error  plots. 

5.4.  The  Peaceman-Rachford  ADI  method  for  parabolic  problems. 

As  a  final  example  we  will  derive  intermediate  boundary  conditions  for  the  Peaceman- 
Rachford  method  (1.42)  by  viewing  this  as  a  time-split  method  for  the  problem  = 
+  «w  with  the  splitting  (1.43),  We  consider  the  problem  on  the  unit  square  0  <  z  < 
1,  0  <  y  <  1  and  assume  that  Dirichlct  boundary  data  is  supplied  at  all  points  on  the 
boundary.  Since  U*  is  differenced  only  in  the  z-direction  in  (1.42),  we  need  to  generate 
intermediate  boundary  data  only  on  the  boundaries  z  =  0  and  z  —  1.  We  will  consider 
only  the  boundary  at  z  =  0.  The  other  boundary  is  completely  analogous. 

The  given  boundary  data  is 


u(0,  y,  t)  =  g[y,  t).  (5.27) 

We  seek  to  determine  U*0m  u*(0 ,mh,tn  +  k)  where  u*  is  the  solution  to  the  first 
subproblem  from  (1.43), 


ut  si11**  +  Vyy)  +  ^k(uSi>„  uvvw)' 

with  inital  conditions 

«*(*.  V,  *n)  =  «(*,  V,  <»),  V  z,  y. 


(5.28) 

(5.29) 
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Differentiating  (5.28)  shows  that 

=  Kvs«s  "h  ^u**yv  "i"  uyyyy)  ^ 

and  so,  proceeding  as  in  Chapter  4, 

u*(0,mh,tn  +  *)  =  u*(0  ,mh,tn)  +  ku*t  (0,  mh,  tn)  +  $*3u«t(0  ,mh,tn)  +  0(k3) 

=  «*  +  +  «^)  +  “  «Cyyy)l 

+  i*2!^1*****  +  2u«yy  +  uyyyy)  +  °(k)\  +  0(k3). 

In  view  of  the  initial  conditions  (5.20),  we  can  replace  u*  by  u  everywhere  on  the 
righthand  side  of  the  last  equation.  By  expanding  u(0,  mh,  tn  +  k/2)  about  u(0 ,mh,t„) 
and  comparing  this  with  what  results  above,  we  find  that 


u*(0,mMn  +  *)  =  t»(0  ,mh,tn  +  k/2)  +  Jfc2(u„,i(0,  mh,tn) 

—  «vvyv(0,  mh,  tn))  +  0(fc3). 


(5.30' 


We  retain  0(Jfc3)  global  accuracy  provided  the  boundary  conditions  have  0(k2)  local 
accuracy.  We  can  thus  neglect  the  0(k3)  terms  in  (5.30)  and  take 

U'0>m  =  «( 0,  mh,  tn  +  k/2)  = 


where  pJJ,+I/3  =  g{mh,tn  +  k/2).  Higher  order  accuracy  at  the  boundary  can  be  obtained 
by  including  the  next  term  as  well.  The  u,,,*  term  cannot  be  calculated  directly  from 
the  given  boundary  conditions  (5.27).  However,  from  the  original  differential  equation 
we  find  that 

—  Cfxxxx  4"  2tf xxyy  +  ®yyy yi 
Utyy  —  Uxxyy  +  u*yyy> 

and  so 

««  2tt|yy  =  USxxx  “  Uyyyy  • 

Since  ti«  and  u(vy  can  both  be  computed  along  the  boundary,  we  can  use  this  expression 
in  place  of  urXx*  —  uyyyy  (3-30),  giving  the  0(Jfc3)  boundary  conditions 


u0,m  =  simh>  ln  +  fc/2)  +  \k3{glt(mh,  tn)  -  2gtvv[mh,  tn)).  (5.31) 


It  is  interesting  to  note  that  if  we  approximate  the  derivatives  of  g  appearing  in  (5.31) 
by  0{k3)  accurate  finite  differences  of  gJJ,  and  at  the  meshpoints,  we  obtain 

V'0,m  =  +  iC")  ~  \kD+vD^(gZ+l  -  »") 

=  Kl  -  i*D+vD_y)C+‘  +  &1  +  i*/?+yD_y)C- 

These  arc  precisely  the  standard  boundary  conditions  for  the  Peaceman-Rachford  method 
as  derived  by  Fairweather  and  Mitchell(19|  using  different  techniques. 

For  irregular  regions  it  is  not  in  general  possible  to  replace  the  z-  and  y-derivatives 
in  (5.30)  by  tangential  and  time  derivatives  which  can  be  evaluated  directly  from  the 
given  boundary  data.  However,  the  expansion  (5.30)  is  still  valid  and  one-sided  finite 
differences  could  be  used  to  approximate  the  fourth  derivatives. 
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